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The  question  of  the  stability  of  steady  state  solutions  in 
geophysical  fluid  flows  is  addressed  through  qualitative  analysis 
and  quantitative  examples.  The  inviscid  linear  stability  theory  of 
stratified  shear  flows  and  the  solution  of  the  stability  problem 
using  normal  modes  and  Fourier-Laplace  transforms  are  discussed. 

Two  numerical  examples  are  used  to  illustrate  the  relationship 
of  various  physical  parameters  to  the  stability  of  the  system  and 
to  trace  the  development  of  the  instability  for  short,  Intermediate 
and  long  times.  The  examples  are  (1)  a  two  layer  fluid  of  infinite 
extent  with  application  to  the  air-sea  interface  and  (2)  a  two- 
layer  fluid  having  a  free  surface  and  finite  depth  with  application 
to  a  salt  wedge  estuary. 

The  initial-value  problem  is  solved  using  a  power  series 
expansion  for  short  times,  superposition  of  modes  for  intermediate 
times  and  asymptotic  analysis  for  long  times.  The  asymptotic 
expansion  applicable  in  non-conservative  systems  is  compared  with 
the  approximate  solution  using  ray  techniques,  which  are  valid  in 
conservative  systems,  and  analytic  continuation  of  the  eigenvalues 
into  the  complex  wavenumber  plane. 


I 


r 


> 


I 


¥ 


> 


TABLE  OF  CONTENTS 


LIST  OF  TABLES 
LIST  OF  FIGURES 
ACKNOWLEDGEl-IENTS 

1.  Introduction 

2.  Initial-Value  Problems  in  Stratified  Shear  Flows 

Equations  of  Motion 
Boundary  and  Initial  Conditions 
Solution  of  the  Problem 
The  Initial-Value  Problem 

3.  Two-Layer  Problem 

Description  of  Problem 
Initial-Value  Problem 
Initial  Period 
Asymptotic  Analysis 

4.  Salt-Wedge  Estuary 

Description  of  Problem 
Mathematical  Description 
Eigenvalues 

Finite  Time  Series  Analysis 
Asymptotic  Analysis 

5.  Summary 
REFERENCES 

Appendix  A.  Fourier-Laplace  Transform  Solution  of  the 
Initial-Value  Problem 


PAGE 

iv 

V 

lx 

1 

6 

6 

11 

15 

21 

25 

25 

34 

35 
39 
41 
41 
44 
53 
57 
59 
65 
67 

70 


» 


Appendix  B.  Asymptotic  Analysis 
Appendix  C.  Tables  and  Figures 


LIST  OF  TABLES 


Page 

Table  I.  Parameters  for  Unbounded  Two-Layer  Problem  83 

Table  II.  Parameters  for  Figures  13  to  22.  89 

Table  III.  Parameters  for  Salt  Wedge  Estuary  Profiles.  100 

Figures  23  to  33. 

Table  IV.  Notation  131 


tv 


LIST  OF  FIGURES 


Page 


Figure 

1. 

Schematic  of  unbounded  two-layer  example. 

26 

Figure 

2. 

Schematic  of  a  longitudinal  section  of  a 
salt-wedge  estuary. 

42 

Figure 

3. 

Schematic  of  downstream  velocity  profile 
in  a  salt-wedge. 

42 

Figure 

4. 

Schematic  of  linearized  profiles  for  salt- 
wedge  example. 

43 

Figure 

5. 

Stability  function  (-Q)  vs.  wavenumber  for 
data  set  la. 

84 

Figure 

6. 

Velocity  of  stability  boundary  vs.  wave- 
number  for  data  set  la. 

84 

Figure 

7. 

Stability  function  (-Q)  vs.  wavenumber  for 
data  set  Ib. 

85 

Figure 

8. 

Stability  function  (-Q)  vs.  small  wave- 
number  for  data  set  Ib. 

85 

Figure 

9. 

Real  part  of  frequency  vs.  wavenumber  for 
data  set  Ic. 

86 

Figure 

10. 

Imaginary  part  of  frequency  vs.  wavenumber 
for  data  set  Ic. 

86 

Figure 

11. 

Initial  distortion  of  Gaussian  pulse  of 
standard  deviation  X  =  0.1. 

87 

Figure 

12. 

Asymptotic  expansion  of  Gaussian  pulse  of 
standard  deviation  X  =  0.1. 

88 

Figure 

13. 

Eigenvalues  for  data  set  A.  No  shear. 

90 

Figure 

14. 

Eigenvalues  for  data  set  B.  No 
stratification 

91 

Figure 

15. 

Eigenvalues  for  data  set  C.  Combined 
profiles. 

92 

Figure 

16. 

Eigenvalues  for  data  set  D.  Top  boundary 
far  from  interface. 

93 

V 


Figure 

17. 

Eigenvalues  for  data  set  E.  Bottom  boundary 
far  from  interface. 

94 

Figure 

18. 

Eigenvalues  for  data  set  F.  Positive 
surface  shear. 

95 

Figure 

19. 

Eigenvalues  for  data  set  G.  Negative 
surface  shear. 

96 

Figure 

20. 

Eigenvalues  for  data  set  H.  Interface  in 
lower  half  of  profile. 

97 

Figure 

21. 

Eigenvalues  for  data  set  I.  Low  discharge. 

98 

Figure 

22. 

Eigenvalues  for  data  set  J.  High  discharge. 

99 

Figure 

23. 

Eigenvalues  for  Profile  1,  Ap  =  0.002. 

101 

Figure 

24. 

Eigenvalues  for  Profile  1,  Ap  =  0.02. 

102 

Figure 

25. 

Eigenvalues  for  Profile  1,  Ap  =  0.2. 

103 

Figure 

26. 

Eigenvalues  for  Profile  2,  Ap  =  0.0. 

104 

Figure 

27. 

Eigenvalues  for  Profile  2,  Ap  =  0.02. 

105 

Figure 

28. 

Eigenvalues  for  Profile  2,  Ap  =  0.2. 

106 

Figure 

29. 

Eigenvalues  for  Profile  3,  Ap  =  0.0. 

107 

Figure 

30. 

Eigenvalues  for  Profile  3,  Ap  =  0.02. 

108 

Figure 

31. 

Eigenvalues  for  Profile  3,  Ap  =  0.2. 

109 

Figure 

32. 

Eigenvalues  for  Profile  4,  Ap  =  0.02. 

110 

Figure 

33. 

Eigenvalues  for  Profile  5,  Ap  =  0.02. 

111 

Figure 

34. 

Time  series.  Initial  exponential  depth 
decay.  Profile  1.  a  =  0.5.  Six  normal 
modes . 

112 

Figure 

35. 

Time  series.  Initial  exponential  depth 
decay.  Profile  1.  ci  =  0.5.  (2N  +  6)  modes. 

112 

Figure 

36. 

Time  series.  Initial  exponential  depth 
decay.  Profile  2.  5  =  0.5.  (2N  +  6)  modes. 

113 

Vi 


Figure 

37. 

Time  series.  Initial  Gaussian  centered  at 
density  interface.  Profile  2.  5=  0.5. 

(2N  +  6)  modes. 

113 

Figure 

38. 

Time  series.  Initial  exponential  depth 
decay.  Profile  3.  a  =  0.5.  {2N  +  6)  modes. 

114 

Figure 

39. 

Line  of  saddle  points.  Most  growing  mode. 
Profile  1. 

115 

Figure 

40. 

Contours  of  3wj^/3u)j.  in  complex  wavenumber 
plane.  Profile  1. 

115. 

Figure 

41. 

Contours  of  frequency  in  complex  wavenumber 
plane.  Profile  1. 

-16 

Figure 

42. 

Wavenumber  vs.  x/t  along  line  of  saddle 
points.  Profile  1. 

117 

Figure 

43. 

Frequency  vs.  x/t  along  line  of  saddle 
points.  Profile  1. 

118 

Figure 

44. 

F  =  (a^x/t  -  0)^)  vs.  x/t  along  line  of  saddle 
points.  Profile  1. 

119 

Figure 

45. 

Asymptotic  expansion.  Initial  flat  spectrum. 
Profile  1,  z  =  0.0,  t  =  10.0. 

120 

Figure 

46. 

Asymptotic  expansion.  Initial  flat  spectrum. 
Profile  1,  z  =  -0.1,  t  =  10.0. 

120 

Figure 

47. 

Asymptotic  expansion.  Initial  flat  spectrum. 
Profile  1,  z  =  -0.2,  t  =  10.0. 

121 

Figure 

48. 

Asymptotic  expansion.  Initial  flat  spectrum. 
Profile  1,  z  =  -0.5,  t  =  10.0. 

121 

Figure 

49. 

Asymptotic  expansion.  Initial  flat  spectrum. 
Profile  1,  z  =  0.0,  t  =  100.0. 

122 

Figure 

50. 

Asymptotic  expansion.  Initial  flat  spectrum. 
Profile  1,  z  =  -0.1,  t  =  100.0. 

122 

Figure 

51. 

Asymptotic  expansion.  Initial  flat  spectrum. 
Profile  1,  z  =  -0.2,  t  =  100.0. 

123 

Figure 

52. 

Asymptotic  expansion.  Initial  flat  spectrum. 
Profile  1,  z  =  -0.5,  5  =  100.0. 

123 

vil 


» 

Figure  53. 

Asymptotic  expansion.  Initial  flat  spectrum. 
Second  mode.  Profile  1,  z  =  0.0,  t  =  10.0. 

124 

Figure  54. 

Asymptotic  expansion.  Initial  flat  spectrum. 
Second  mode.  Profile  1,  z  =  -0.5,  t  =  10.0. 

124 

r 

Figure  55. 

Asymptotic  expansion.  Initial  flat  spectrum. 
Second  mode.  Profile  1,  z  =  0.0,  t  =  100.0. 

125 

Figure  56. 

Asymptotic  expansion.  Initial  flat  spectrum. 
Second  mode.  Profile  1,  z  =  -0.5,  t  =  100.0. 

125 

f 

Figure  57. 

Asymptotic  expansion.  Initial  Gaussian 
spectrum.  Profile  1,  z  =  0.0,  t  =  10.0, 

A  =  0.1. 

126 

1 

Figure  58. 

Asymptotic  expansion.  Initial  Gaussian 
spectrum.  Profile  1,  z  =  0.0,  t  =  100.0, 

A  =  0.1. 

126 

Figure  59. 

Ray  mathematics  expansion.  Initial  flat 
spectrum.  Profile  1,  z  =  -0.1,  t  =  10.0. 

127 

¥ 

Figure  60. 

Ray  mathematics  expansion.  Initial  flat 
spectrum.  Profile  1,  z  =  -0.2,  t  =  10.0. 

127 

Figure  61. 

Ray  mathematics  expansion.  Initial  flat 
spectrum.  Profile  1,  z  =  -0.5,  t  =  10.0. 

128 

» 

Figure  62. 

Ray  mathematics  expansion.  Initial  flat 
spectrum.  Profile  1,  z  =  -0.1,  t  =  100.0. 

129 

Figure  63. 

Ray  mathematics  expansion.  Initial  flat 
spectrum.  Profile  1,  z  =  -0.2,  t  =  100.0. 

29 

» 

Figure  64. 

Ray  mathematics  expansion.  Initial  flat 
spectrum.  Profile  1,  z  =  -0.5,  t  =  100.0. 

130 

Acknowledgements 


1  would  like  to  express  my  sincere  appreciation  to  Professor 
William  0.  Criminale,  Jr.  for  the  continued  guidance,  support  and 
encouragement  in  bringing  this  work  to  successful  completion. 

I  am  also  particularly  grateful  to  my  friends  and  fellow 
students  who  have  offered  suggestions,  advice  and  encouragement 
throughout  my  studies. 

This  research  was  supported  by  the  Air  Force  Office  of 
Scientific  Research  under  Grant  AFSOR  74-2579. 


Introduction 


In  the  theoretical  study  of  geophysical  fluid  flows  an  attempt 
is  rade  to  find  a  steady-state  solution  to  the  problem.  The  steady- 
state  solution  may  be  the  result  of  solving  the  governing 
equations  expressing  the  balance  of  forces  in  the  physical  problem 
or  :t  may  be  the  result  of  considering  time  averaged  motion.  In 
either  case  it  is  of  considerable  interest  to  know  if  the  steady- 
state  solution  persists  or  if  the  flow  is  dynamically  unstal'le  and 
evolves  into  a  time  dependent  flow. 

The  technique  employed  to  answer  this  question  is  the  linear 
theory  of  the  hydrodynamic  stability  of  parallel  flows.  The 
problem  is  posed  as  a  combined  boundary-value  and  initial-value 
problem.  A  time  dependent  perturbation  expansion  is  made  about 
the  steady-state  solution.  The  perturbations  are  assumed  to  be 
small  compared  to  the  mean  flow  and  the  equations  are  linearized. 
Traditionally  the  emphasis  has  been  on  the  solution  of  the  boundary- 
value  or  eigenvalue  problem.  The  eigenvalues  are  evaluated  to 
determine  if  the  system  has  an  exponentially  growing  mode.  If  an 
exponentially  growing  mode  is  found,  then  the  system  is  unstable 
to  small  perturbations  and  the  steady-state  solution  does  not 
persist . 

Completing  the  solution  by  solving  the  initial-value  problem 
is  more  difficult  because  it  requires  solving  for  the  eigen¬ 
functions  as  well  as  the  eigenvalues.  For  some  problems  it  is 
possible  to  determine  the  eigenvalues  without  being  able  to  solve 
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for  the  eigenfunctions.  However,  the  solution  to  the  initial 
value  problem  is  of  Interest  because  it  not  only  indicates  that  the 
solution  is  time  dependent,  but  also  determines  the  structure  of 
the  time  dependence. 

For  most  geophysical  fluid  flows  it  is  possible  to  formulate 
a  very  general  set  of  equations  describing  the  balance  of  forces 
and  the  applicable  boundary  conditions  and  initial  conditions. 

The  resulting  set  of  equations  forms  a  well-posed  mathematical 
problem,  but  is  analytically  insoluble.  To  overcome  this  diffi¬ 
culty  approximations  are  made  at  the  expense  of  limiting  the  appli¬ 
cability  of  the  resulting  solution. 

Justification  for  these  approximations  and  the  resulting 
limitations  can  be  found  in  the  literature.  For  example,  the 
evaluation  of  the  eigenvalues  of  a  three-dimensional  problem  has 
been  reduced  to  the  solution  of  an  equivalent  two-dimensional 
problem  through  a  transformation  due  to  Squire  (1933),  however, 
the  full  three-dimensional  structure  of  the  initial-value  problem 
is  sacrificed. 

Simplification  of  the  systems  to  be  considered  is  the  result 
of  the  work  of  many  investigators.  The  use  of  the  linearized  Euler 
equations  as  an  approximation  to  the  linearized  Navier-Stokes 
equations  for  flows  at  high  Reynolds  number  and  small  wavenumbers 
is  demonstrated  in  work  by  Case  (1961)  and  Lin  (1967).  This  is  an 
Important  simplification  because  it  reduces  the  order  of  the 
differential  equation  to  be  solved.  Drazin  (1958,  1961)  has  shown 


that  it  Is  possible  to  approximate  continuous  density  and  velocity 
profiles  by  discontinuous  or  piecewise  continuous  profiles. 

The  eigenfunctions  of  the  boundary  value  problem  form  a  set  of 
normal  modes  which  are  used  as  a  basis  for  expanding  the  initial 
conditions.  In  general,  however,  the  eigenfunctions  do  not  form  a 
complete  set  and  it  is  not  possible  to  expand  arbitrary  initial 
conditions  as  a  sum  of  the  eigenfunctions.  Case  (1960)  and 
Pedlosky  (1963)  demonstrated  that  in  addition  to  the  discrete 
spectrum  of  normal  modes  there  are  also  continuous  spectra  of 
modes  arising  from  the  singular  solutions  of  the  Inviscid  differen¬ 
tial  equations.  Evaluation  of  the  initial-value  problem  must 
include  these  solutions  to  be  valid. 

Following  this  it  was  pointed  out  by  Banks,  Drazln  and 
Zaturska  (1976)  that  the  modes  of  a  parallel  inviscid  stratified 
fluid  fall  into  five  classes,  some  of  which  may  be  absent  for  a 
given  flow.  These  five  classes  are: 

(1)  a  finite  set  of  non-singular  unstable  modes 

(2)  a  conjugate  finite  set  of  non-singular  damped 
stable  modes 

(3)  a  finite  set  of  singular  stable  modes,  each  having 
a  branch  point  and  being  the  limit  of  an  unstable 
mode 


(4)  a  discrete  set  of  non-singular  stable  modes,  and 
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(5)  a  continuous  set  of  singular  stable  modes. 

The  complex  nature  of  stratified  shear  flows  is  indicated  by 
the  findings  of  several  investigators.  It  is  generally  agreed 
that  increasing  the  stable  stratification  and  the  presence  of 
boundaries  stabilize  the  flow.  However,  Maslowe  and  Kelly  (1971) 
showed  that  under  some  conditions  increasing  the  stable  strati¬ 
fication  leads  to  destabilization  of  the  flow.  Hazel  (1972)  and 
Lalas,  Einaudi  and  Fua  (1976)  have  demonstrated  the  destabiliz¬ 
ing  effects  of  boundaries.  Davey  (1977)  showed  that  increasing 
the  curvature  of  the  mean  velocity  profile  increases  the  range  of 
unstable  modes. 

Two  different  mathematical  models  are  used  to  demonstrate 
features  of  the  stability  and  initial-value  problems.  The  choice 
of  these  examples  is  dictated  by  the  existence  of  analytic 
solutions  to  demonstrate  features  of  the  initial-value  problem 
while  maintaining  relevance  to  geophysical  flows.  The  first 
system  is  an  infinite  two-layer  flow  with  uniform  velocity  and 
density  in  each  layer  and  surface  tension  at  the  interface.  The 
eigenvalues  for  this  system  exhibit  a  limited  range  of  unstable 
modes  with  a  definite  maximum  growth  rate.  This  system  is  used 
to  demonstrate  the  three-dimensional  development  of  an  initial 
Gaussian  pulse  for  short  tim.  and  in  the  asymptotic  limit.  Results 
from  this  example  have  application  to  disturbances  of  the  air-sea 


interface. 


The  second  system  is  modeled  on  features  of  the  salt-wedge 
estuary.  The  fluid  is  of  finite  vertical  extent;  bounded  above 
by  a  free  surface  and  below  by  a  solid  wall.  The  density 
structure  is  two  layers  and  the  continuous  velocity  profile  is 
approximated  by  three  piecewise-linear  segments.  This  model  is 
used  to  investigate  the  effects  of  stratification,  shear  and 
boundaries  on  stability.  The  difference  between  a  correct 
asymptotic  exapnsion  of  the  initial-value  problem  of  a  non¬ 
conservative  system  and  the  approximate  ray-mathematics  expansion 


is  demonstrated. 


2. 


Initial-Value  Problems  in  Stratif ied-Shear  Flows 


Equations  of  Motion 

The  initial-value  problem  corresponds  to  a  description  of  the 
motion  of  the  fluid  derived  from  the  physical  laws  governing  the 
motion  and  the  constraints  on  that  motion  in  the  form  of  boundary 
conditions  and  initial  conditions.  The  governing  equations  for 
inviscid  and  incompressible  fluid  flow  are: 

Conservation  of  momentum  - 

+  u  •  Vu  =  —  Vd  -  VG  (2.1) 

ot  P  — 

Conservation  of  mass  - 

1^  +  u  .  Vp  =  0  (2.2) 

Incompressibility  - 

V  •  u  =  0  (2.3) 

The  motion  is  described  as  a  combination  of  a  mean  flow  plus  a 
perturbation  about  the  mean.  The  mean  quantities  of  the  flow  are 
time  independent  and  for  the  purpose  of  this  analysis  are  only  a 
function  of  the  vertical  spatial  dimension.  The  perturbation 
quantities  are  functions  of  time  and  the  three  spatial  dimensions. 
It  is  assumed  that  the  mean  motion  is  entirely  horizontal. 
Utilizing  these  restrictions  the  variables  in  the  equations  of 
motion  take  the  following  form: 
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-► 

u  = 

(u*. 

v' 

.  w*) 

u’  = 

IN 

+ 

u(x,  y,  z, 

t) 

y'  = 

Viz) 

+ 

v(x,  y,  z. 

t) 

w'  = 

w(x, 

z,  t) 

e'  = 

p(z) 

P(x,  z, 

t) 

Po 

T 

P 

o 

= 

0(z) 

+ 

e(x,  y,  z. 

t) 

P’  = 


1 


°  p(z’)dz'  +  p(x,  y,  z,  t) 


(2. A) 


The  equations  of  motion  for  the  small  perturbations  on  the  mean 
flow  are  obtained  by  substituting  the  functional  form  of  the 
variables  indicated  in  eq.  (2.4)  into  the  governing  equations 
(2.1,  2.2  and  2.3).  The  terms  in  the  equations  are  collected  by 
order;  the  zeroth  order  terms  govern  the  mean  motion.  These  terms 
are  subtracted  to  obtain  the  equations  for  the  higher  order  terms. 
The  remaining  equations  are  linearized  in  terms  of  the  perturbation 
variables.  These  are: 
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2 

V  p  = 


'dz  3x  dz  -  -  3z 


(2.8) 


Another  equation  Involving  the  Laplacian  of  the  pressure  can 
be  derived  by  taking  the  Laplacian  of  the  vertical  component  of  the 
momentum  equation. 


"  3l 


a  a  a  2  dy  .  dl^  3w 

ro  .r,  d  „r-d  .  -  d'v 


-  f—  +  Tj  —  +  V  V  w  -  2  ( —  —  +  —  — 1 
^3t  -  3x  -  3y  -  Mz  3x  dz  3y 


3y'  3z 


d  y  „  d  V  .  2 

( —  +  — -  T— )  w  +  g7  e 

dz2  dz2  -  - 


(2.9) 


Elimination  of  the  pressure  between  the  vertical  derivative  of  eq. 
(2.8)  and  eq.  (2.9)  yields 


a  a  a  2  d^Ud^V 

(I7  +  y  I-  +  |-)v  w  -  (—  f-  +  —  I-)  w 

^3t  -  3x  -  -  ^  2  3^  j_2  3r  - 


-gyj  e 


(2.10) 


g2  32 

Where  =  —-2  +  —3 

h  3x'^  3y‘‘ 


To  eliminate  the  density  from  eq.  (2.10)  first  take  the  horizontal 
Laplacian  of  the  equation  for  the  conservation  of  mass,  eq.  (2.6) 
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(2.13) 


dO  o2 

-  5  31  'h  ^ 


0 


Boundary  and  Initial  Conditions 

The  boundary  conditions  to  be  applied  in  solving  the  differen¬ 
tial  equation  are  determined  by  the  physical  problem  under  con¬ 
sideration.  For  the  salt  wedge  estuary  problem  it  is  assumed  that 
the  horizontal  extent  is  much  greater  than  the  vertical  extent. 
Because  of  the  great  difference  in  the  horizontal  to  vertical 
length  scales,  the  side  boundaries  have  little  effect  on  the 
motion  at  the  scales  under  investigation.  Therefore,  the  x  and 
coordinates  are  considered  to  be  infinite  in  extent.  The  approp¬ 
riate  boundary  conditions  at  x  — >■  +»  and  ^  — >  ±“>  are  that  w  and 
the  first  derivatives  of  w  in  the  horizontal  directions  remain 
finite. 

The  boundary  conditions  in  the  vertical  direction  are  the 
free  surface  condition  at  the  surface  and  a  no  flux  condition  at 
the  bottom. 

The  free  surface  is  a  constant  pressure  surface  and  the 


boundary  condition  is 
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0P_ 

Dt 


0 


(2.14) 


This  condition  can  be  expressed  in  terms  of  the  vertical  velocity. 
The  pressure  is  expressed  as  the  sum  of  three  components;  the 
atmospheric  pressure,  th  hydrostatic  pressure  and  the  perturba¬ 
tion  pressure.  The  functi  lal  dependence  of  each  of  these 
components  is  indicated  in  the  following  equation. 


E'  (^» 

“  “a  ”  “ 

+  g[0(z)  +  e(x,  z,  t)]  •  [n(x,  y,  t)  -  z]  (2.15) 


+  g(x,  y,  z,  t) 


where  n(x,y,t)  -  z  =  0  is  the  equation  of  the  free  surface. 
Expanding  the  free  surface  condition  yields 


g(0  +  9)  w 


+  8(0  +  0)  (|j  +  ^  +  £|-)  „ 
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+  B(n  -£)%  +  %•»  E|^)  0 


^  90  ,  N  J.  /■  .90 

■•■  S  ^  S^D  "  5^  i! 


(2.16) 


+  u 


^  ^  -vlA 


9x  -  9^ 


;  P  +  w 


9p 
-  9z 


=  0 


Let  the  free  surface  be  defined  as 


+ 


gOn) 


=  0 


(2.17) 


This  is  the  zeroth  order  surface  boundary  condition  and  is  sub¬ 
tracted  from  eq.  (2.16).  The  remaining  terms  are  linearized  and 
rewritten  utilizing' conservation  of  mass  and  the  definition  of  a 
material  surface  at  z  =  q.  The  linearized  surface  boundary 
condition  becomes 


(2.18) 


To  the  same  order  of  approximation  this  condition  is  applied  at 
z  =0,  rather  than  at  the  perturbed  surface  z  =  n. 

The  required  form  of  the  boundary  condition  is  an  equation  in 
w  only.  To  eliminate  the  pressure  from  the  above  equation,  use  the 
horizontal  momentum  equations  and  continuity.  This  is  accomplished 
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by  operating  on  the  x-component  of  the  momentum  equations  by 


3x 


and  the  y-component  by 


3_ 


Add  the  resulting  equations  and  solve  for 


V —]  72 
-  3y^  h 


This  is  substituted  into  eq .  (2.18)  after  taking  the  horizontal 


Laplacian.  The  resulting  boundary  condition  is 


9w 

Jz 


a  a 

-  '•dz  3x 


dV  3 
dz  3^^ 


w 


g  0  V2  w  =  0 


(2.19) 


at  z 


0. 
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The  boundary  condition  at  the  solid  bottom  is  simply  no  flow 
through  the  bottom  or 


=  0 


(2.20) 


at  z  =  -H ,  the  depth  of  the  bottom  boundary. 


The  initial  conditions  for  the  problem  are  the  specification 

of  the  values  of  w(x,y,z,0)  and  (x.^.z.O),  or  some  combination 
3w 

of  w  and  at  time  t  =  0.  Such  conditions  will  be  considered  in 
greater  detail  in  conjunction  with  specific  applications. 

Solution  of  the  Problem 

Two  methods  of  solution  for  the  differential  equation  speci¬ 
fied  in  eq.  (2.13)  are  to  be  considered,  namely,  normal  modes  and 
Fourier-Laplace  transforms. 

The  first  method  considered  is  the  normal  modes  solution. 

This  is  essentially  a  separation  of  variables  method  which  assumes 
that  the  horizontal  space  and  time  dependence  of  the  solution  is 
wave-like  and  has  the  form: 


-  “  w(z)[Ae^°-  +  Be^°-]e^^*^-  ^ 


(2.21) 


Substitution  of  this  assumed  form  for  the  solution  into  the  partial 
differential  equation  (2.13)  results  in  an  ordinary  second  order 
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differential  equation  for  w(z) .  This  equation  is 

2  a2 

(kU  +  ZV  -  a)  (^2  -  -  £2)  v7 

—  —  .  dz  — 

d20  d^V  . 

-  (k^  +  ZV  -  o)  (k-^2  +  w  (2.22) 

dO 

-  g  ^  (k^  +  Z^)w  =  0 

Now  consider  the  Fourier-Laplace  transform  method  of 
solution.  This  technique  utilizes  a  double  Fourier  transform  in 
X  and  y  and  a  Laplace  transform  of  the  partial  differential 
equation  (2 . 13) . 

The  double  Fourier  transform  is  defined  by  the  transform  pair 


F(w)  =  w(k,£,z,t)  =  ^2^72"  /”  /”  ^  ^^^dxdy 

'  ^  _ m  _ an 


(2.23) 


F-^w)  =  w(x.y,z,t)  =  r  r  w(k.^z,t)e^(’"2  + 


dkd*, 


The  Laplace  transform  of  Fourier-transformed  function  is 


defined  by  the  transform  pair 
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L(w) 


=  w(k,!!,,z,s)  =  j  w  (k.H.z.t 


t)e  dt 


(2.24) 


L  ^(w)  =  w(k,^,z,t)  = 


e+1® 


St 


w(k,ll,z,s)e  -  ds 


’e-i“ 


Hence,  the  function  w(k,J.,z,s)  is  defined  by  the  triple  integral 


FL(w) 


w(k,il,z,s) 


(2.25) 


(2tt)‘ 


/OO  ^  f<30 

Jn  J-oo  J_a 


w(x,y,z,t)e 


-i(kx  +  iy) 


e  dxdydt 


The  known  properties  of  Laplace  and  Fourier  transforms  of  deriva¬ 
tives  which  are  necessary  for  the  transform  of  eq.  (2.13)  are: 


9'^w\ 


3x 


(-ik)*^  F(w) 


'■at'' 


=  s  L 


(w)  -  w  1^, 


0+ 


L  (w)  -  s  w 


aw 


-  't=o+  at  '  t=o+ 
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The  boundary  conditions,  eqs.  (2.19)  and  (2.20)  are  also  Fourler- 
Laplace  transformed  and  take  the  form; 

dw 

(s  +  ikU  +  ilV)^ 

dj/  dF  ^ 

-  (s  +  iky  +  lAT)  i  (2.28) 

+  gO  (k^  +  Jl^)w=0  atz  =  0 

and 

w  =  0  at  z  =  -H  (2.29) 

It  can  be  seen  that  the  Fourier  Laplace  transformation  of  the  partial 
differential  equation  in  w  leads  to  an  inhomogeneous  second  order, 
ordinary  differential  equation  in  the  transformed  variable 
w  (k, i,z,s).  Making  the  substitution  s  =  -io  in  eq.  (2.27)  and 
comparing  the  homogeneous  part  with  eq.  (2.22)  it  can  be  seen  that 
the  derived  governing  equation  for  the  two  methods  is  of  the  same 
form. 

Since  the  differential  equation  is  of  the  same  form  for  the  two 
methods,  that  given  in  eq.  (2.22)  will  be  used  to  formulate  the  non- 
dimensional  problem.  The  equation  will  be  non-dimensionalized  on  a 
velocity  scale,  u  ,  a  length  scale,  h  ,  and  a  density  scale,  p  . 


The  non-dimensional  variables  are  defined: 
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(“.TT) 

(x.y.z) 

(u,v,w) 

(U,V) 

P 

u> 

0 

0 


(kh  ,  dh  ) 
o  o 

(x/h  ,  y/h  ,  z/h  ) 

—  o  *•  o  —  o 

(u/u  ,  v/u  ,  w/u  ) 

—  o  —  O  "■  o 

(i//u  ,  v/u  ) 

“0-0 

p/u^ 

-  o 

o  h  /u 
o  o 

P/Po 

p/p^ 


Making  the  substitutions  into  eq.  (2.13)  produces  the  non- 
dimensional  form  of  the  governing  differential  equation: 


(aU  +  yV  -  li)} 


-  iaU  +  -  w)  (at/"  +  yV”)  w 


(2.30) 


where 


and 


+  J  (o^  +  Y^)  w  =  0 


N^ho^ 

- ^  =  Richardson  No. 

u  ^ 
o 


=  -g0  '  =  Brunt-Vdisala  frequency. 


Utilizing  the  Squire  transformation  it  can  be  shown  that  the  solu¬ 
tion  to  the  three-dimensional  problem  posed  above  is  equivalent  to 
the  solution  of  a  two-dimensional  problem  with  appropriately  scaled 


parameters. 
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The  Squire  transformation  is  defined  by  the  following  two 
definitions. 


=  a  + 

all  =  all  +  yV 


(2.31) 


Using  these  definitions  in  eq.  (2.30)  the  equivalent  two-dimensional 
problem  is  given  by 


{all  -  w)^(^^j  -  a^)  w  -  a(aU  -  u)  U"  w 

+  Jw  =  0  (2.32) 

where 

J  = 

Applying  the  Squire  transformation  to  the  boundary  conditions,  eq. 
(2.19)  and  eq.  (2.20),  completes  the  specification  of  the  two- 
dimensional  problem.  The  boundary  conditions  for  the  equivalent 
two-dimensional  problem  are 


(aU  -  uj)2  w’  -  [a2g0  +  aU'  (ai/  -  m)]  w  =  0  (2.33) 

at  z  =  0 

and 

w  =  0  at  z  =  -H/h 

o 

The  Initial-Value  Problem 

Solution  of  the  initial-value  problem  using  the  normal  modes 


» 
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approach  requires  solution  of  the  eigenvalue  problem  specified  by 
the  differential  equation  for  w(z),  eq,  (2. 22),  subject  to  the 
boundary  conditions  given  by  eqs.  (2.19)  and  (2.20).  The  eigen¬ 
functions  can  be  written  as 


N  . 

w(a,Y,z,t)  =  I  w^(a,Y,z) [A^e 
n=l 


-10)  t  .  „  io)  t] 
n  +  B  e  n  ^ 
n 


(2.34) 


where  the  n  specifies  the  mode  and  corresponds  to  the  eigenvalue 
given  by 


0)  =0)  (a.Y,"!) 

n  n 


(2.35) 


To  apply  the  initial  conditions  to  the  problem,  first  take  the 

Fourier  transform  or  find  the  Fourier  series  of  the  initial  vertical 

velocity  and  acceleration.  Let  w  (x,y,z,0)  and  —  w  (x,y,z,0)  be 

O  d  t  O 

the  initial  conditions  and  define  the  Fourier  transforms  of  these 
conditions  as 


w^(a,Y,z,0) 


■^w^(a,Y»z,0) 


w  (x,y,z,0)  e  ^^^dxdy 

o 

(2.36) 

^  (x.y.z.O)  g-i(ax  + 


Now  expand  the  transforms  of  the  initial  conditions  in  terms  of  the 
eigenfunctions  of  the  normal  modes  solution,  i.e.. 
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N  > 

w  (a,Y,z,0)  =  (a,Y.z)[A  +  B  ] 

o  n  n  n 

(2.37) 

3  '  N  - 

—  w^(a.Y.z.O)  =  I  -  la)^w^(a.Y.z)[A^  -  B^] 
n=l. 


The  coefficients  A  and  B  are  determined  from  eqs.  (2.37)  and  the 
n  n 

final  solution  to  the  initial-value  problem  is  found  by  taking  the 
inverse  Fourier  transform  of  the  normal  modes  solution  which  is  of 
the  form 


w(x,y,z, t) 


=  f  [  I  w^(a,Y,z)  [A^e 
^  ^  n=l 


(2.38) 


+  B  e^“ntj 
n 


i(ax  +  Yy) 


dadY 


The  solution  to  the  initial-value  problem  using  the  Fourier-Laplace 
transform  method  begins  by  solving  for  the  Green's  function  solution 
of  the  differential  equation  given  by  eq.  (2.27)  which  is  rewritten 
as 

•2 

(s  +  iaU  +  -  y^)  w 

-  (s  +  iaU  +  ±yV)  (laU"  +  i.yV")  w  (2.39) 

«  A 

+  g0'  (a^  +  Y^)  w  =  P/(a,Y.z;s) 


where  the  function 
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f/(a,Y.z;s) 


y2) 


3w 

at 


o"^ 


+  I 


s  +  21(aU  +  yV)]  [ 


dz2 


(2.40) 


-  i(ai/"  +  yV^’)  w 


specifies  the  initial  conditions.  The  functions  w  and  dw/dt  are  the 
Fourier  transforms  of  the  vertical  velocity  and  acceleration  as 
specified  in  eqs.  (2.36).  The  Green's  function  solution  to  eq. 
(2.39)  is  discussed  in  Appendix  A  and  is  of  the  form 


G(z,Zp)  =  G(z,z^;s,a,Y) 


(2.41) 


The  solution  to  the  initial-value  problem  in  terras  of  the  Green's 
function  is  specified  formally  by  the  inverse  Fourier-Laplace 
transform: 


“  “>  e+i®  o 


w(x,y,z,l:)  .^1  j  j  I 

—CO  —00  ^00— 

St  i(ax  +  Yy) 

•  e  e  ■' 


G(z,z^)  V(a,Y,z^;s)  dz^ 


(2.42) 


dsdadY 


» 


The  correspondence  between  this  solution  and  the  normal  modes 
solution  is  also  discussed  in  Appendix  A. 


3. 


Two-Layer  Problem 
Description  of  Problem 

The  first  example  of  an  initial-value  problem  is  a  two  layer 
system  consisting  of  two  fluids  of  different  density  in  uniform 
motion  relative  to  each  other  with  surface  tension  at  the  interface. 
It  is  assumed  that  the  fluids  are  inviscid  and  incompressible.  The 
illustration  (Fig.  1)  shows  the  general  features  of  the  system. 

The  density  of  the  bottom  layer,  p^,  is  greater  than  the  density  of 
the  upper  layer,  p^.  The  system  is  statically  stable.  The  analysis 
assumes  that  the  velocity  of  the  upper  layer  is  +U  and  that  of  the 
bottom  layer  is  -U.  The  essential  feature  is  the  velocity  dif¬ 
ference  W.  A  Galilean  transformation  could  be  made  to  place  the 
zero  of  the  velocity  at  any  appropriate  value  without  changing  the 
following  analysis. 

The  scaling  parameters  for  this  problem  are  based  on  the 
magnitude  of  the  velocity  Z/,  the  gravitational  constant  g  and  the 
density  of  pure  water,  1  gm/cm^.  The  velocity  scale  u^  =  £/  and  the 
length  scale  h^  =  Z/^/g.  The  non-dimensional  variables  are  as 
defined  in  Section  2.  The  non-dimensional  coefficient  of  surface 
tension,  x,  is  defined  as: 

T  =  T  u^h  . 

-  o  o 

In  the  preceding  section  the  differential  equation  governing 
the  motion  of  an  inviscid  incompressible  stratified,  shear  flow  is 
given  by  eq.  (2.8) 
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Figure  1,  Schematic  of  unbounded 
two-layer  example. 
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(all  +  ~  -  Y^)  w 

dz 

-  (aU  +  yV  -  ui)  (aU"  +  yV")  w 

+  J  (ci2  +  Y^)  w  =  0 


(3.1) 


For  this  example  the  differential  equation  can  be  greatly 
simplified  -  the  horizontal  velocity  in  the  y-direction  is  zero 
and  the  velocity  in  the  x-direction  is  constant  in  each  of  the 
layers.  Also,  the  density  is  constant  in  each  of  the  two  layers. 
This  means  that  the  derivatives  of  U  with  respect  to  z  are  zero 
and  J  =  0.  The  differential  equation  for  w  in  each  of  the  layers 
is  simply 


dz2 


(a^  +  Y^)  w  =  0 


(3.2) 


or 


dz2 


2 

a^w 


0 


which  has  as  solutions 

/  \  A  OtZ  ,  _  ”CtZ 

w(z)  *  Ae  +  Be 


The  solution  in  each  region  will  be  designated  by  a  subscript  1  or  2. 

It  is  assumed  that  the  fluid  is  of  infinite  extent  and  the 
boundary  conditions  are  applied  at  z  ±  “  rather  than  at  z  =  0  and 
z  =  -h.  The  Interface  between  the  two  fluids  at  z  ■*  0  is  an 
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internal  boundary  which  divides  the  solution  into  two  regions. 
Matching  conditions  to  be  applied  to  the  solutions  in  the  two 
regions  must  be  derived. 

The  two  conditions  which  are  to  be  applied  are  the  kinematic 
and  dynamic  boundary  conditions.  The  kinematic  condition  is  that 
the  interface  is  a  material  surface  specified  by 


z  = 


.  i(ax  +  vy  -  mt) 
n  =  d  e  ^ 


(3.4) 


and  is  subject  to  the  conditions  that 


Dt 


(3.5) 


The  linearized  form  of  the  kinematic  condition  is 


•^  +  1/-^  =  w 
3t  3x  + 


at  z  =  0 


(3.6) 


and 


in  -  yla  =  w 

3t  3x  - 


at  z  =  0 


where  0  and  0  indicates  the  value  at  0  approached  from  above  and 
below  respectively. 

The  dynamic  condition  is  that  the  difference  in  pressure  across 
the  interface  is  sustained  by  surface  tension. 

Brenoulli's  equation  for  the  pressure  can  be  written  as 

p=-|^-gz--|u*  u  (3.7) 
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where  is  the  velocity  potential  of  the  flow  defined  by 

u  =  -V<p  (3.8) 

The  boundary  condition  is 


The  continuity  equation  gives  the  differential  equation  which  the 
velocity  potential  must  satisfy. 

V^<f,  =  0  (3.12) 
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which  gives 
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^^(z)  =  b^e““^ 
i’jCz)  =  32®°*^ 


(3.17) 


Substituting  the  functional  form  of  ij)  and  n  into  the  conditions 

(a),  (b)  and  (c)  gives  the  following  set  of  linear  algebraic 
equations  to  be  solved  for  the  eigenvalues  ,  u). 


(a)  i(a)  -  al/)d  +  ab^  =  0 

(b)  i(a)  +  ay)d  -  aa^  =  0  (3.18) 

(c)  lg(0,  -  0,)  -  xa^ld  +  i0  (u  -  ay)b  -  i0  (u  +  aU)a  =  0 

12  1  12  2 

The  determinant  of  coefficients  of  these  three  equations  is  the 
eigenvalue  relation  and  can  be  written  as  the  solution  of  a  quad¬ 
ratic  equation: 


-g  cosij)  VLQ 

(01  +  ©2) 


-4a^  cos^</< 


(01  +  ©2)2 


ggA© 


3 

g^T 


(©1  +  ©2)  (©1  +  ©2) 


1/2 


(3.19) 


where 


A©  =  ©2  -  ©1 
g  =  g  cos  Ip 
Y  =  g  sin  Ip 


tan  4)  =  y/g 


The  stability  boundary  for  this  problem  is  obtained  from  con¬ 
sidering  the  argument  of  the  radical  in  eq.  (3.19).  Let  this  be 
represented  by  Q, 


-  cos^tjy  0  ) 

Q  =  a  - _ - _ U.  +  _ m _ 

^  1(01  +  ©2)  (©1  +  ©2^^  (©1  +  ©2^J 


(3.20) 


If  Q  is  negative,  w  is  complex  and  the  flow  can  be  unstable,  but  if 
Q  is  positive  (u  is  always  real.  The  stability  boundary  is  deter¬ 
mined  by  the  locus  Q  =  0 .  The  range  of  instability  can  be  determined 
by  solving  the  quadratic  equation  in  brackets. 

(©^  +  ©2)  (0j  +  0^)2  (©j  +  ©2) 

If  the  roots  of-  eq.  (3.21)  are  complex,  the  flow  is  stable  or 
unstable  everywhere  depending  on  the  sign  of  Q.  However,  if  the 
roots  are  real,  they  specify  the  points  at  which  Q  changes  sign  and 
thus  the  regions  of  instability.  Figures  5,  7  and  8  show  the  values 
of  (-Q)  for  two  different  cases. 

The  values  of  the  parameters  used  in  evaluating  Q  for  Fig.  5 

are: 

T  •=  74.  cm^/sec2 
Pj  =  1.293  X  10  2  gm/cm^ 
pj  “  1.02  gm/cm^ 

U  ■  100.  cm/sec 
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The  parameters  for  this  case  correspond  to  an  air-water  interface. 
The  roots  of  eq.  (3.21)  for  this  case  are  complex  and  the  value  of  Q 
is  positive.  The  two-layer  problem  with  this  set  of  parameters  is 
stable  for  all  wave-numbers.  The  threshold  velocity  for  instability 
is  324.405  cm/sec  at  a  wavenumber  of  3.673/cm.  The  value  of  the 
threshold  velocity  as  a  function  of  wavenumber  is  shown  in  Fig.  6. 
The  wavenumber  range  for  instability  is  given  by  the  points  of 
intersection  of  a  line  of  constant  velocity  and  the  threshold 
velocity  curve. 

The  second  case  corresponds  more  closely  to  an  internal  fluid 
boundary.  The  surface  tension  is  small  and  the  density  difference 
between  the  two  fluids  is  small.  The  parameters  are 

T  =  1.0  cm^/sec^ 
pj  =  1.0  gm/cm^ 

P2  =  1.02  gm/cm^ 
y  =  10.  cm/sec 

Figure  7  is  a  graph  of  the  value  of  (-Q)  as  a  function  of  wave- 
number  and  indicates  the  regions  of  instability  for  waves  at  dif¬ 
ferent  angles,  i/>,  to  the  mean  velocity.  Figure  8  is  the  same 
.  function,  but  indicates  the  stable  (-Q  <  0)  region  of  the  flow  at 
small  wavenumbers  (large  wave  lengths). 

The  eigenvalues  for  this  set  of  parameters  are  shown  in  Figs.  9 
and  10.  Figure  9  shows  the  real  part  of  the  eigenvalue  u.  It  can 
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be  seen  that  the  real  part  of  the  frequency  is  nearly  constant  over 
the  range  of  instability.  The  imaginary  part  of  w  is  shown  in 
Fig.  10.  As  in  all  inviscid  problems,  the  complex  eigenvalues  occur 
in  complex  pairs.  For  every  growing  mode  there  is  also  a  damped 
mode.  The  range  of  instability  is  clearly  shown  in  this  graph. 

Initial-Value  Problem 

The  solution  to  the  initial-value  problem  is  given  by  (2.38)  as 


w(x,y,z,t)  =  [  [  I  w^(a,Y,z)  [A^e 


(3.22) 


with  the  initial  conditions 


N 


w  (a.Y.z.O)  =  w  (a.Y.z)  [A  +  B  1 
o  n  n  n 

n=l 


and 


3w  N  . 

(a,Y,z,0)  =  J  w  (a,Y.z)  [-iw  A  +  im  B  ] 
,  n  n  n  n  n 

n=l 


3t 


(3.23) 


As  an  example  of  an  initial-value  problem,  the  flow  described  by 
the  second  set  of  parameters  in  the  preceding  section  is  used.  The 
Initial  conditions  prescribed  are  chosen  as 


w^(a,Y.z,0)  = 


3w 

(a,Y,z,0)  =  0 


■(2^P~ 


-a2z2 


(3.24) 
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This  is  the  Fourier  transform  of  a  Gaussian  pulse  in  amplitude 
centered  around  the  density  interface.  Substituting  the  prescribed 
initial  conditions  into  eq.  (3.23)  determines  the  value  of  and 
Solving  and  making  the  substitution  into  eq.  (3.22)  gives  the 
result  that 


w(x,y,z,t) 


00  00 

.  Del  I  I 

n=l  •'  ‘ 


(a.Y.z.O)  e 


l(ax  +  Yy  -  oJt) 


dady 


(3.25) 


Initial  Period 

To  determine  the  motion  of  the  pulse  for  small  times,  the 
function  is  expanded  in  a  Taylor  series  about  time  t  =  0.  The 
following  technique  follows  that  described  by  Criminale  (1960).  The 
Taylor  series  expansion  takes  the  form 


,  \  V  9^  (  f  '  f  \  l(c»x+Yy-  wt) 

w(x,y,z,t)  =  I  ^ j  w^(a,Y.z)e 

•dadYl^O 


.00  .00 


(3.26) 


The  Fourier  Integral  in  eq.  (3.25)  can  be  broken  into  two  Integrals 


and  rewritten  as 
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(3.31) 
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The  evaluation  of  only  one  of  these  integrals  is  necessary 
since  the  other  is  just  the  complex  conjugate. 

The  information  needed  is  the  form  of  the  time  derivative  of  W 
evaluated  at  time  t  =  0. 


t=0 


f  ~  .  „k  i(ax  +  Yy) .  ■ 

w  (-imj  e  ■'  dady 

;  o 

-o 


(3.32) 


If  0)  can  be  approximated  by  a  polynomial  these  integrals  can  be 
evaluated  explicitly.  Each  term  is  of  the  form 


3t 


=  (-io))' 


i(ax  +  Yy)j  j 
we  ■'  dady 

o 


t=0 


(3.33) 


This  form  is  obtained  by  recalling  the  expression  for  the  Fourier 
transform  of  a  derivative 

F(^)  =  (-ia)''F(w) 

3x 

The  linear  operator  w  in  eq.  (3.33)  has  the  same  functional  form  as 
(i)(a,y)  with  a  replaced  by  (-1  3/3x)  and  y  by  ("i  3/3y). 

Approximate  the  eigenvalues  u  by  the  following  polynomials : 


(I)  =  Re(a))  =0)  +  a,o  +  a,y^ 

r  ro  1  2 ' 

(3.34) 

=  Im(a))  =  -  bj  (a  - 

Also  note  that  the  integral  in  eq.  (3.33)  is  half  the  initial 


» 


condition  in  real  space 


M 
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2  w^Cx.y.z) 


w  oo 

-  f  f  ~  ^  (“X  Yy) .  3 

-  J  J  ”o  ®  ^^''dadY 


(3.35) 


Define  the  complex  linear  operator  L  such  that 


L  =  -1(4)  =  L  +  IL 
-  r  1 


32  .  .  02 


"  “i  [w  +  2a  b ,  -r—  -  a  — —  1 
1  ro  o  1  ^ 


3y2 


(3.36) 


3x  2  3y2 

Substituting  into  the  Taylor  series  eq.  (3.31)  takes  the  form 


w(x,y,z,t)  =  j  —  [l^  +  L 


k=0 


(3.37) 


Expanding  this  series  and  keeping  only  real  terms  the  initial 
disturbance  takes  the  form 


w(x,y,z,t)  =  w  +  t  (w  ) 
o  r  o 


+  t2[L2  -  l2] 


(3.38) 


Expanding  the  operator  L  and  keeping  only  the  linear  term,  the 
expansion  valid  for  small  time  is  given  by 
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w(x,y,2,t)  =  {1  +  t[(w  -  b  a2) 

lO  i  o 


b.x'^ 


(3.39) 


+  a 


1  x/A^  -  bj/A^  +  — - b^/A^ 


bjy^ 

A** 


-]} 


(x^  +  y2)/x2 
4TrA2 


The  initial  spreading  and  distortion  of  the  initial  conditions 
given  by  eq.  (3.24)  are  shown  in  Fig.  11.  The  initial  Gaussian 
pulse  has  a  standard  deviation  equal  to  0.1.  Fig.  11a  shows  the 
initial  pulse  at  t  =  0.  The  distorted  pulse  is  shown  at  t  =  0.05 
and  t  =  0.75  in  Fig.  11b  and  Fig.  11c,  respectively. 

Asymptotic  Analysis 

The  form  of  the  asymptotic  expansion  of  the  solution  to  the 
initial-value  problem,  eq.  (3.25),  is  derived  in  Appendix  B.  This 
expansion  is  given  by  eq.  (B.19) 


w(x,y,2,t) 


'V  Re 


,  .  i(a*x  +  y.y  -  uj.t) 

w^(a*,Y*,z,o)  e  '  *  *  ' 


3^u) 

[32(0 

2 

3a2 

’W 

3a3y 

1/2 


(3.40) 


where  and  y*  ^re  solutions  of  the  simultaneous  equations; 
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(3.41) 

y/t  -  |^(a*.Y*)  =  0 

where  to,  a  and  y  ai'e  considered  to  be  complex  variables. 

Evaluation  of  the  expansion  for  the  same  Gaussian  pulse  with 
standard  deviation  X  =  0.1  is  shown  in  Fig.  12.  The  expansion  in 
Fig.  12a  is  at  t  =  20  and  that  in  Fig.  12b  at  t  =  60. 

The  continued  spreading  and  distortion  of  the  pulse  and  the 
appearance  of  dominant  wavelengths  are  apparent  in  the  series  of 
figures  11  and  12  showing  the  development  of  the  pulse  from  t  =  0 


to  t  =  60. 


Description  of  Problem 

A  salt  wedge  estuary  is  used  as  the  second  physical  example  in 
this  analysis. 

A  salt  wedge  estuary  is  a  highly  stratified  estuary  consisting 
of  nearly  homogeneous  upper  and  lower  layers  separated  by  a  strong 
halocline.  The  brackish  upper  layer  is  primarily  river  runoff.  The 
"salt  wedge"  is  the  layer  below  the  density  interface  consisting 
of  a  saltwater  intrusion  from  the  ocean.  Salt  wedges  occur  in 
shallow  estuaries  with  high  river  discharge  and  relative] v  weak 
tidal  currents. 

The  general  profile  of  a  salt  wedge  is  depicted  in  Fig.  2  and 
a  general  velocity  profile  in  Fig.  3.  The  velocity  profile  is  such 
that  there  is  no  net  flow  in  a  stationary  salt  wedge.  Three  regions 
are  indicated  on  the  schematic  of  the  wedge  profile.  These  are 
defined  by  the  depth  of  the  halocline.  The  outer  region  has  the 
density  interface  near  the  surface.  The  inner  regime  is  defined  by 
a  halocline  near  the  bottom  with  the  central  region  having  a  halo¬ 
cline  near  the  center  of  the  water  column. 

The  velocity  profile  is  nearly  two-dimensional  and  is  treated 
as  such  in  the  following  analysis.  In  addition,  the  scale  of 
perturbations  to  be  considered  are  such  that  the  x-dependencc 
(downstream)  of  the  velocity  and  density  structure  is  considered 
only  in  terms  of  the  inner,  central  and  outer  regime.  At  each 
cross-section  the  velocity  and  density  profile  will  be  similar  to 
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Figure  2.  Schematic  of  a  longitudinal  section 
of  a  salt-wedge  estuary. 
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Figure  3.  Schematic  of  downstream  velocity 
profile  in  a  salt-wedge. 
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Figure  4.  Schematic  of  linearized  profiles 
for  salt-wedge  example. 
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those  shown  in  Fig.  4.  The  parameters  h^  and  are  scaled  on  the 
total  depth  of  the  fluid  and  the  gravitational  constant  g. 

h  =  H 
o 

“o" 

Mathematical  Descviption 

The  density  structure  is  modeled  as  a  two-layer  system  with 
density  pj  in  the  upper  layer  of  depth  d2»  The  lower  layer  is  of 
density  P2  of  thickness  h-d2.  The  smooth  velocity  profile  is 
modeled  by  a  piecewise  linear  profile  consisting  of  three  segments. 
The  first  segment  is  given  by 

Uj^  =  Uj  +  Uj  z  ;  0  ^  z  ^  -d^ 

where 


Uj  =  =  velocity  at  the  surface 

dGi 

u'  =  —  =  shear  in  layer  I 

1  dz 


u  -  u 
s  _ m 


Similarly,  the  middle  segment  of  the  velocity  profile  is  given  by 


‘II 


'II 


+  u 


II 


-dj  ^  z  i  -dj 


where 
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“ll  “  (d3  -  d^) 

.  ~  ^ 
“ll  “  (dj  -  dj) 


and  the  bottom  segment  is  given  by 


‘III 


=  “ill  “in  " 


-dj  ^  z  ^  -h 


where 

“ill  ~  (h  -  dg) 

.  “b 

“in  (h  -  dg) 

A  restriction  on  the  velocity  profile  is  that  it  goes  to  zero  at  the 
bottom. 

These  approximations  divide  the  vertical  extent  of  the  model 
into  four  regions.  In  each  region  the  mean  properties  of  the  flow 
are  described  by  constant  density  and  constant  shear.  There  are 
three  internal  boundaries  defined  by  discontinuities  in  the  density 
or  velocity. 

The  differential  equation  governing  the  flow  and  the  external 
boundary  conditions  are  given  by  eqs.  (2.30),  (2.31)  and  (2.32). 
Repeating  the  equations  and  using  the  definition  of  J  gives 


iaU  -  0))^  ^  ~  ^  ~  80*  aw  =  0  (^*1) 


dw  \  a^gO 


icxU  -  ta)^  {aU  -  w) 


w  =  0  at  z  ®  0 


(4.2) 


w  =  0  at  z  *  -h 


(A. 3) 


The  differential  equation  simplifies  Immediately  on  substi¬ 
tuting  the  assumed  profiles  of  mean  velocity  and  density  to 


(aU  -  w  =  0 


The  solutions  of  this  equation  are  subject  to  the  boundary  condi¬ 
tions  eq.  (4.2)  and  eq.  (4.3)  and  to  matching  conditions  at  the 
discontinuities  of  the  mean  profiles. 

The  applicable  matching  conditions  are  derived  by  integrating 
the  complete  differential  equation  given  by  eq.  (4.1)  across  the 
discontinuity.  The  matching  condition  at  depth  z  =  -d  is  given  by 
the  following  equation. 


(aU  -  10)2 


A 

d^w  ,  '2 f 

^dz-a2J 


(aU  -  u))2  wdz 


-a-i-e 

■  f  ' 

-  a  w 


(aU  -  (o)  y"dz  -  ag  wG'dz  =0 


(4.4) 


Evaluation  of  the  integrals  at  a  given  depth  will  be  slmpll- 

A 

fied  by  requiring  that  the  vertical  perturbation  velocity,  w(z),  be 
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continuous  and  noting  that  the  mean  velocity  Is  a  continuous  func¬ 
tion  of  depth,  but  the  shear  Is  discontinuous. 

Evaluating  eq.  (4.4)  at  z  =  d^,  proceeds  as  follows. 
Evaluating  the  Integrals  and  taking  the  limit  of  small  e  produces 


^  *"(i  * 

(“i’m  -  ‘  -  O 


-d 


-  otw  (-djXayjn  -  w)[^]  '  -0 


=  0 


-d 


1 


which  can  be  solved  for  dw/dz  to  give 


(4.5) 


aw(-d,) 

(— ]  =  — : -  (u'  -  u'  )  at  z  =  -d,  (4.6) 


Evaluating  the  Integrals  In  eq.  (4.4)  at  the  other  discontinuities 
leads  to  the  remaining  two  matching  conditions. 

At  z  =  -d2 


-d. 


-ag(02  -  0j)  w(-d2) 

(a  Up  -  u)^ 


(4.7) 


where 


UP  =  Ujj  -  ujj  d2 


and  at  z  =  -d. 
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>^1 


a  w(-dj) 

s  -I  — 

-dj  -  w) 


(4.8) 


The  eigenvalue  problem  for  the  salt  wedge  estuary  is  given  by 
the  solution  to  the  following  set  of  equations: 


Wj  =  0 


(4.9) 


which  has  solutions  of  the  form 


Wj(z)  =  aj  sinh  (az)  +  ^ j  cosh  (az) 


(4.10) 


where  i  indicates  the  layer  in  which  the  solution  is  valid. 

A 

The  conditions  on  the  solutions  Wj(z)  are: 


dw,  '  u| 

'•>  TT  ■  >«(/  -  (5h'  -  ■»)'  °  at  ^  -  0 


(b)  Wj  =  w 


at  z  =  -d. 


(c) 


dwj  dw2 
dz  dz 


a  w 


(ay  -  u) 


m 


at  z  =  -d. 


(d)  w,  =  w. 


at  z  =  -d. 


dw- 


(e) 


dWg  oig(®2  "  ®i)  ^2 


(4.11) 


dz  dz 


(a  Up-  ui) 


at  z  =  -d. 


(f)  W3  =  w^  at  z  =  -dj 
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(g) 


dw^ 

dz 


aw. 


(ay,  -  ui)  ^'^li  ■ 


at  z  =  -dj 


(h)  =  0  at  z  =  -h 

The  result  of  substituting  the  solution  Wj  into  eqs.  (4.11)  is  a  set 
of  six  linear  equations  in  the  unknown  coefficients  aj  and  b j . 

The  determinant  of  coefficients  of  this  set  of  linear  equations  is 
the  dispersion  relation  for  the  system  and  can  be  written  as  a  sixth 
degree  polynomial  in  the  frequency,  w. 

+  4*310^  +  <J>2U^  +  tjU)  +  =  0  (4.12) 

where  the  <tj  are  complicated  transcendental  functions  of  a,  the 
mean  velocity,  shear  and  density. 

The  eigenfunctions  take  the  form: 

w^(z)  =  a^  (r  sinh  az  +  a  cosh  az) 

W2(z)  =  a^  sinh  az  +  b2  cosh  az 

(4.13) 

Wj(z)  =  Sj  sinh  az  +  b^  cosh  az 
w^(z)  =  a^  sinh  a(z  +  h) 

where 

a2  =  aj[r  +  \|)  cosh  adj(r  sinh  adj  -  a  cosh  ad^)] 
b2  =  aj[a  +  i/i  sinh  ad  (r  sinh  ad^  -  a  cosh  ad^)] 
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“^2  “^2^  ~  *’2'^  cosh^adj 

^3  “  ^2^  sinh^ad^  +  b2[l  -  ((>  sinh  ad  cosh  adj] 


a,.  = 


[-33  slnh  adg  +  bg  cosh  ad^J/sinh  a(h  -  d^) 


r  =  TT 


aS0 


au' 


<1'  = 


(aU^  -  a))2  (a[/^  -  a,) 

(u^  -  up 


(fii/  -  0)) 
m 


ag(G2  -  Oj) 
(ai/p  -  a))2 


The  next  step  In  the  solution  of  the  problem  posed  in  eq.  (4.4) 
is  to  evaluate  the  solution  for  the  values  of  (o  =  aU.  Because  the 
coefficient  (aU  -  u)  is  squared  there  are  two  separate  solutions 
possible  for  each  u  =  aU.  These  solutions  form  the  continuous 
spectra  of  modes  which  Case  (1960)  showed  are  necessary  to  solve 
the  initial  value  problem.  They  are  given  by  the  solutions  to  the 
following  differential  equations.  The  first  is 


d^w 

s 

dz2 


6(z  +  z^) 


(4.14) 


where  6(z  +  z  )  is  the  Dirac  delta  function  and  -z  is  the  depth  at 
o  o 

which  the  velocity  U  is  equal  to  w/a. 


The  second  is 
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» 


dz2 


w  =  {'(z  +  z  ) 
b  o 


(4.15) 


where  6'(z  +  z  )  Is  the  derivative  of  the  Dirac  delta  function  with 
o 

respect  to  z.  Making  the  substitution  Wj^  =  the  equation  to  be 

solved  becomes 


iia  _ 

dz2 


a^q 


6(z  +  z  ) 
o 


(4.16) 


The  solution  of  these  two  equations  proceed  in  the  same  manner  as 
solution  of  eq.  (4.9)  except  that  there  is  an  additional  condition 
to  be  applied  at  z  =  ~z^. 

A 

Integrating  the  equation  for  w^,  eq.  (4.14)  from  z  =  -z^  -  e 
to  z  =  -z^  +  e  gives  the  following  condition  to  be  applied  at 


z  =  -z 


-z  +e  o”' 
f  o  d^v; 


-z  +e 


-z  +e 


dz2 


dz 


-  I  w^  dz  =  I  6(z  +  z^)  dz  (4.17) 


-z  -e 
o 


-z  -e 
o 


-z  -e 
o 


Taking  the  limit  as  e  -*•  0  and  requiring  the  function  w  to  be 
continuous  at  z  =  -z^  gives 


Z“-Z 


+ 

o 


1 


(4.18) 


Similarly,  the  condition  to  be  applied  to  q(z)  at  -z^  is 
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1 


(4.19) 


where  the  requirement  is  that  q  be  continuous  at  However, 

from  the  definition  w^^  =  ^  it  can  be  seen  that  the  jump  condition 
eq.  (4.19)  is  actually  a  jump  in  given  by 

w  (-Z  '*’)  -  w.  (-Z  ")  =  1 
bo  bo 


The  solution  of  the  homogeneous  equation  (4.4)  for  the  salt 
wedge  is  divided  into  four  regions  with  boundary  conditions  at  the 
surface  (4.11a)  and  the  bottom  (4.11h)  and  matching  and  jump  con- 

A  A 

ditions  given  by  eq.  (4.11b-g).  The  solution  w  and  w  to  the 

&  D 

two  auxiliary  equations  (4.14)  and  (4.15)  must  also  satisfy  the 
conditions  given  by  eqs.  (4.11a-h).  However,  the  region  is  divided 

A  A 

into  five  layers  for  the  solutions  of  w^(z)  and  w^(z).  The  position 
of  the  additional  internal  boundary  depends  on  the  depth  z  =  -z^ 
which  is  defined  by  U(-z^)  -  m/a  and  can  occur  anywhere  in  the 
region  0  ^  z  ^  -h.  The  solutions  have  the  form 


w 

ai 


w 


34 


=  A[r  sinh  oz  +  a  cosh  az] 

=  32  sinh  az  +  b2  cosh  az 

=  33  sinh  az  +  b^  cosh  az 

=  34  sinh  az  +  b4  cosh  az 

=  3  sinh  a(h  +  z) 

5 


(4.21) 
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and 


dq. 

(Z)  . 


where 


qj  =  A' [a  slnh  az  +  P  cosh  az] 

I  -  1  - 

q^  =  cosh  az  +  sinh  az 

1  -  I  - 

=  a^  cosh  az  +  sinh  az  (4.22) 

I  ~  I  ~ 

=  a^  cosh  az  +  b^^  sinh  az 

I 

q^  =  a^  cosh  a(h  +  z) 

The  relationships  among  the  coefficients  A,  a^  and  b^  and 

I  I 

A',  a^  and  b^  depend  on  the  depth  -z^  relative  to  “dj,  -a^  and  -d^ 

and  are  different  for  -z  in  each  layer. 

o 

Eigenvalues 

m,  A  A 

The  value  of  to  (a)  for  the  auxiliary  functions  w  (z)  and  w,  (z) 
are  given  by  a)(a)  =  ai/(z^)  and  form  a  continuous  spectrum  of  neutral 
modes.  These  correspond  to  the  continuous  spectrum  of  neutral 
modes  which  arise  in  the  solution  of  the  initial-value  problem 
using  Laplace  transforms.  In  the  normal  mode  solutions  of  the 
initial-value  problem  these  provide  additional  modes  for  fitting 
arbitrary  initial  conditions  (cf.  Chimonas,  1977). 
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The  eigenvalue  problem  specified  by  the  differential  eq.  (4.9) 
and  the  conditions  (4.11a-h)  lead  to  the  sixth  degree  polynomial 
given  by  eq.  (4.12), 

^  '  i 

F(a,io)  =  ^  (a)  0)  (4.23) 

j=0  ^ 

itj  (a)  =  transcendental  functions  of  a. 

In  general  F(a,w)  is  a  complex  function  of  two  complex 
variables  a  and  w.  However,  the  stability  problem  is  usually 
solved  as  a  temporal  problem  by  considering  a  real  or  as  a  spatial 
problem  in  which  u  is  real.  Solution  of  the  temporal  problem 
results  in  u  being  given  as  a  function  of  a  as  in 

u)  =  f(a)  (4.24) 

For  this  specific  example  the  temporal  problem  results  in  a  dis¬ 
crete  spectrum  of  six  eigenmodes  which  are  the  solutions  of  the 
polynomial  given  in  eq.  (4.12). 

The  spatial  problem  has  an  eigenvalue  relation  of  the  form 

a  =  g(o))  (4.25) 

which  for  this  example  is  the  solution  of  a  non-linear  transcen¬ 
dental  equation.  As  a  result  of  this  nothing  can  be  said  in 
general  about  the  number  of  modes,  except  that  the  spectrum  is 


discrete. 
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The  temporal  eigenvalues  of  the  salt-wedge  problem  are  investi¬ 
gated  to  determine  the  effects  of  stratification,  surface  shear,  top 
and  bottom  boundaries  and  transport  of  the  water  column.  The 
numerical  parameters  for  this  study  are  given  in  Table  II. 

The  eigenvalues  for  parameter  set  A  shown  in  Fig.  13  shows  the 
stable  modes  associated  with  the  two-layer  density  profile  with  no 
mean  motion.  The  unstable  mode  in  the  homogeneous  fluid  with  the 
velocity  profile  of  set  B  is  shown  in  Fig.  14.  This  profile  has  no 
shear  at  the  surface.  The  eigenvalues  for  parameter  set  C,  which 
is  a  combination  of  A  and  B,  shows  the  interaction  of  the  two  pro¬ 
files  resulting  in  four  unstable  modes  (Fig.  15). 

The  effect  of  the  top  boundary  on  the  eigenvalues  is  shown  in 
set  D  which  moves  the  top  boundary  far  from  the  density  interface. 
This  suppresses  two  modes  as  shown  in  Fig.  16.  Reversing  the 
procedure  and  moving  the  bottom  boundary  far  from  the  density 
interface  is  given  by  set  E  (Fig.  17).  Comparing  this  case  with 
set  C  (Fig.  15)  shows  that  one  mode  is  suppressed  by  removing  the 
effect  of  the  bottom  boundary.  A  small  positive  shear  at  the 
surface  (Set  F,  Fig.  18)  alters  the  relative  position  of  the  four 
unstable  modes  and  tends  to  decrease  the  magnitude  of  the  two  modes 
associated  with  the  surface.  The  effect  of  a  positive  shear  at  the 
surface  is  to  decrease  the  curvature  of  the  velocity  profile.  Data 
Set  G  (Fig.  19)  represents  a  negative  shear  at  the  surface.  This 
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effect  is  to  increase  the  curvature  of  the  velocity  profile  and 
Increase  the  amplitude  of  the  two  unst'’  le  surface  modes. 

Moving  the  density  interface  to  a  position  near  the  bottom  as 
in  Set  H  (Fig.  20)  has  the  effect  of  suppressing  all  but  one  mode. 

Comparing  Data  Sets  I  and  J  (Figs.  21  and  22)  shows  that 
increasing  the  transport  and  consequently  increasing  the  shear 
causes  an  increase  in  the  value  of  the  unstable  mode  and  intro¬ 
duces  an  additional  unstable  mode. 

The  parameters  of  the  profiles  representing  salt-wedge 
estuaries  are  given  in  Table  III.  Profiles  1,  2  and  3  represent 
the  outer  region  with  no  surface  shear,  positive  surface  shear  and 
negative  surface  shear,  respectively.  Profile  4  represents  the 
central  regime  and  Number  5  the  inner  regime. 

The  series  of  Figs.  14,  23,  24  and  25  show  the  effect  on  the 
eigenvalues  of  increasing  the  density  difference  between  the  layers. 
For  the  homogeneous  fluid  there  is  one  unstable  mode.  Two  unstable 
modes  are  added  with  the  two- layer  structure.  As  the  density 
difference  becomes  greater,  the  strength  of  the  instability  increases 
as  well  as  the  frequency.  Figures  23,  24  and  25  show  that  the 
scale  length  of  the  instability  decreases  as  the  density  difference 
increases. 

Profile  2  which  has  a  positive  surface  shear  exhibits  only 
one  unstable  mode  over  the  range  of  density  differences  evaluated. 
Figures  26,  27  and  28  show  that  with  the  decreased  curvature  of  the 
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profile,  the  increased  density  appears  to  stabilize  the  profile  and 
decrease  the  scale  of  instabilities. 

The  negative  shear  of  Profile  3  not  only  increases  the  curva¬ 
ture,  but  also  effectively  introduces  a  second  inflection  point. 
This  profile  has  two  unstable  modes  for  the  homogeneous  fluid  as 
shown  by  the  eigenvalues  in  Fig.  29.  Introducing  a  density  dif¬ 
ference  adds  one  unstable  mode  at  the  scales  under  consideration 
as  shown  in  Fig.  30.  Further,  increasing  the  density  difference 
either  suppresses  that  unstable  mode  or  moves  it  to  a  much  smaller 
scale  as  shown  in  Fig.  31.  Profiles  4  and  5  each  exhibit  two 
unstable  modes  for  the  scales  being  considered  as  shown  in  Figs.  32 
and  33,  respectively. 


Finite  Time  Series  Analysis 

The  solution  to  the  initial-value  problem  is  given  by  eq. 
(2.36)  or  for  the  two-dimensional  problem. 


w(x,z,t) 


00 

r  N  .  - 

=  I 

n=l 


.  -iwnt  ,  _  Imnt,  iax,  ,, 

z)  [A  e  ^  +  B  e  1  e  da  (4.46) 
n  n 


with  initial  conditions  specified  by 


,  -  N  .  > 

w^(a,z)  =  I  w^(a,z)  [A^  + 
n=l 


9w 

if 


N  .  . 

I  -ia)„Wn(a.z)  [A  -  B^] 
n=l 


(4.47) 


r 
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The  initial  conditions  to  be  considered  are  of  the  form 

w^(x,z)  =  f(z)e^“o*  (4.48) 

and 

3w 

if  °  ° 

which  have  transforms 

w  (a,z)  =  f(z)  6(o  -  a  )  (4.49) 

o  o 

and 

3w 

if  =  0 

Substituting  the  initial  conditions  eq.  (4.49)  into  (4.47) 
results  in  a  set  of ‘linear  equations  in  the  coefficients  and  B^. 

N  .  . 

f(z)  =  I  2A  w  (a  ,z)  (4.50) 

-  n  n  o 
n=l 

The  salt-wedge  estuary  problem  has  six  discrete  eigenmodes  plus  two 

continuous  spectra  from  the  singular  solutions.  The  solution  is 

formed  from  the  six  discrete  eigenmodes  and  a  fixed  number,  Nj, 

of  the  continuous  modes.  A  least  squares  fit  at  (2N^  -f  6)  equally 

spaced  positions  z^  is  made  to  eq.  (4.50)  to  determine  the  (2Nj  +  6) 

coefficients  A  . 

n 
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Two  forms  of  the  function  f(z)  are  considered.  One  is  an 
exponential  decaying  from  the  surface  and  the  other  is  a  Gaussian 
pulse  centered  at  the  density  interface.  The  exponential  function 
is  of  the  form 

*/v  2xz/A  . 

f(2)  =  e  ,  X  =  0.1 


and  the  form  of  the  Gaussian  pulse  is 


f(z) 


-(z  +  d2)^/x2 
e 


X  =  0.1 


Figure  34  is  an  expansion  of  the  exponential  at  a  =  0.5  using  only 
the  six  discrete  modes  from  Profile  1.  Including  21  modes  from  each 
of  the  singular  solutions  produces  Fig.  35  which  has  a  difference 
in  phase  and  amplitude  at  a  given  time.  Figures  36  and  37  show  the 
difference  between  the  exponential  and  Gaussian  z-dependence  for 
an  expansion  at  a  =  0.5  using  eigenvalues  from  Profile  2.  At  this 
wavelength  there  are  no  growing  modes  included  in  the  expansion. 

An  expansion  of  the  exponential  at  a  =  0.5  for  Profile  3  is  shown 
in  Fig.  38. 


Asymptotic  Analysis 

The  asymptotic  evaluation  of  the  vertical  velocity  perturbation 
given  by 

00 

-I  f a; c;.z) 

n*l  •* 


w(x,z,t) 


(4.51) 


60 


is  derived  in  Appendix  B.  The  asymptotic  form  of  the  vertical 
velocity  perturbation  is  given  by  eq.  (B.12) 


where  to  (a  )  is  the  most  unstable  eigenmode, 
n 

Figure  39  shows  the  line  of  saddle  points  of  the  most  unstable 

mode  for  an  asymptotic  expansion  using  Profile  1  with  Ap  =  0.02. 

Each  point  on  the  line,  which  is  defined  by  3u)^/9a^  =  0,  represents 

3(0i 

a  different  value  of  x/t.  The  contours  of  t —  near  the  saddle  point 

3otr 

on  the  real  axis  are  shown  in  Fig.  AO.  Figs.  Ala  and  b  show  con¬ 
tours  of  0)^  and  to^  in  the  same  region.  The  values  of  the  wave- 

number  and  frequency  along  the  line  of  saddle  points  are  given  in 
Figs.  A2  and  A3,  respectively.  The  argument  of  the  exponential  in 
the  asymptotic  expansion,  (a^x/t  “<*>„)»  is  plotted  as  a  function 
of  x/t  in  Fig.  AA.  The  region  of  growth  of  the  exponential  is 
clearly  defined  by  the  negative  imaginary  part  of  this  function 
which  extends  approximately  from  -.15  to  .95.  Figs.  A2a  and  A3a 
show  that  over  this  range  the  real  part  of  the  wavenumber  and 
frequency  is  nearly  constant. 

The  asymptotic  expansion  for  a  flat  spectrum  of  wavenumbers  is 
shown  for  t  =  10  at  different  depths  in  Figs.  A5,  A6,  A7  and  A8. 

The  same  expansion  evaluated  at  t  =  100  is  shown  in  Figs.  A9,  50, 
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51  and  52.  At  this  time  the  packet  has  become  peaked  around  the 
most  growing  ray.  For  the  same  initial  conditions  the  expansion 
}  of  the  second  most  growing  mode  is  shown  in  Figs.  53  and  54  at 

t  =  10  and  in  Figs.  55  and  56  at  t  =  100.  At  time  t  =  10,  this 
mode  has  not  fully  developed  into  a  wave  packet,  but  the  wave 
packet  is  obvious  at  t  =  100.  The  amplitude  factors  indicate  that 
this  mode  is  two  orders  of  magnitude  smaller  than  the  dominant  mode 
at  t  =  10  and  over  20  orders  of  magnitude  smaller  at  t  =  100. 

Figs.  57  and  58  show  the  asymptotic  expansion  of  the  solution 
driven  by  a  Gaussian  distribution  of  standard  deviation  0.1.  The 
similarity  of  this  profile  and  that  for  a  flat  spectrum  is  to  be 
expected  because  of  the  nearly  constant  value  of  a  along  the  line  of 
saddle  points. 

The  most  difficult  part  of  the  asymptotic  analysis  is  determi¬ 
nation  of  the  line  of  saddle  points.  The  computation  is  straight¬ 
forward  but  time  consuming  if  u)(a)  is  a  known  analytic  function. 

If  Id  is  known  only  for  real  values  of  a  and  the  first  two  deriva¬ 
tives  of  0)  with  respect  to  a  are  known  or  can  be  approximated, 
simple  ray  mathematics  is  often  used  to  evaluate  the  asymptotic 
form.  However,  Caster  (1978)  warns  that  in  non-conservative  sys¬ 
tems  this  technique  may  lead  to  an  incorrect  result  and  should  be 
used  with  care.  To  demonstrate  the  erroneous  result  which  can 
occur  an  example  is  calculated  to  compare  with  the  correct 


asymptot.'c  expansion. 


I 


Following  the  description  of  the  technique  in  Caster  (1978) 
the  eigenvalue  and  its  first  derivative  are  expanded  about  a  point 
on  the  real  a  axis.  The  wave  number  on  the  real  axis  and  functions 
evaluated  at  that  point  are  designated  by  the  subscript  o. 

The  Taylor  series  expansion  for  the  eigenvalue  u  about  the 
point  on  the  real  axis  a  is  given  by 


.  ,  .du.l,  V  d^o)  . 

u>  =  u)  +(a-a)T^  {o.  -  a  )  -jzrr  +  . . . 

o  o  da  2  o  da2 

o 


(4.53) 


The  expansion  for  the  first  derivative  is 


do)  du  .  ~  x  d'^w  . 

^  =  da  -  “o>  ^ 

o  o 


(4.54) 


The  point  on  the  real  axis  is  such  that 

x/t  =  Re(f  I  ) 
o 

Expanding  eq.  (4.54)  into  real  and  imaginary  parts  gives 


o 


(4.55) 


_  -d(o. 


_  ydoj  X  ,  \  d^bJ  , 

lm(—  )  +  (ct  -  cx^)  ^  +  ... 

o  o 


(4.56) 


If  a  is  to  be  a  saddle  point,  Im(^)  must  equal  zero. 

Equation  (4.56)  can  be  solved  for  the  saddle  point  a. 


f 
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(A. 57) 


It  should  be  noted  that  this  is  valid  only  for  small  values  of 


W.  I 

/  d^o) 

/  da^ 

'  '  o 

/  a 

0 

o 

which  insures  that  the  saddle  point  lies  close  to  the  real  axis. 
The  value  of  the  eigenvalue  at  the  saddle  point  is  evaluated  by 
substituting  eq.  (4.57)  into  eq.  (4.53). 


0).  =  u 
*  o 


(4.58) 


+  ... 


The  form  of  the  asymptotic  expansion  as  given  by  eq.  (B.12)  is 

,it(a^  x/t  -  0)^) 


w(x,z,t)  -v 


w  (“*«z)  e  * 
n  * 


% 

d^u 

1/2 

da^ 

* 

(4.59) 


Making  the  substitutions  from  eq.  (4.55),  (4.57)  and  (4.58),  the 
exponent  becomes 


it (a .  x/t  -  m.)  =  it 


o  x/t  -  u)  +-:r 
o  o  2 


do) . 


dd 


12 


dd2 


(4.60) 
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Substituting  from  eq.  (4.60)  into  (4.59),  the  asymptotic  form, 
assuming  the  validity  of  ray  mathematics,  is 


w(x,z,t) 


it' 

e 

O  x/t-O)  +-^ 
o  o  2 

do)^ 

da 

« 

Vd^to) 

/  da^ 

r 

o 

/ 

o 

d^u) 

1/2 

da2 

o 

(4.61) 


Figures  59-64  show  the  evaluation  of  the  asymptotic  form  given 
by  eq.  (4.61)  for  an  initial  condition  giving  a  flat  spectrum  and 
using  eigenvalues  from  Profile  1.  The  form  of  the  function  indi¬ 
cates  a  caustic.  The  rays  as  defined  above  cross  and  the  disturb¬ 
ance  at  two  or  more  points  map  into  the  same  region. 

Comparing  Figs.  59-64  with  Figs.  45-52  which  are  the  results 
of  the  asymptotic  expansion  for  the  same  initial  conditions  and 
flow  parameters  clearly  indicate  the  care  which  must  be  exercised 
in  using  simple  ray  mathematics  as  an  approximation  to  the 


asymptotic  limit. 


t 


5.  Summary 

The  equations  specifying  the  Inviscid  linear  stability  of 
j  stratified  shear  flows  are  derived  and  the  solutions  using  normal 

modes  and  Fourier-Laplace  transforms  are  discussed.  It  is  shown 
that  solution  of  the  initial  value  problem  by  superposition  of  the 
normal  modes  can  be  accomplished  only  by  including  the  continuous 
spectra  of  singular  solutions  because  the  normal  modes  do  not  form 
a  complete  set.  Examination  of  the  Fourier-Laplace  transform  solu¬ 
tion  indicates  that  it  is  equivalent  to  the  normal  mode  solution 
including  the  singular  solutions. 

The  infinite  two  layer  system  is  used  to  demonstrate  a  power 
series  technique  for  evaluating  the  short  time  distortion  of  an 
initial  pulse.  The  three  dimensional  structure  of  the  asymptotic 
development  of  the  pulse  is  also  shown. 

The  complicated  interaction  of  boundaries,  density  interface, 
shear  and  curvature  of  the  velocity  profile  in  determining  the 

■  stability  of  a  system  is  demonstrated  by  the  eigenvalues  calculated 

i 

I 

for  the  finite  two-layer  piecewise-continuous  velocity  profile 

I  of  Section  A.  Varying  the  parameters  of  the  system  has  led  to 

generalizations  which  may  be  applicable  to  a  salt-wedge  estuary. 

I  Increasing  shear  in  the  velocity  profile  leads  to  more  unstable 

i 

modes.  In  the  physical  system  this  could  occur  during  periods  of 
Increased  river  runoff.  A  downstream  wind  leads  to  a  positive  shear 
In  the  velocity  profile  at  the  surface  and  a  decrease  in  the 
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curvature  of  the  velocity  profile  which  tends  to  stabilize  the 
system.  Conversely,  an  upstream  wind  increases  the  curvature  and 
tends  to  destabilize  the  system.  Over  the  range  of  density  dif¬ 
ferences  included  in  the  calculations,  it  appears  that  increasing 
the  stratification  generally  destabilizes  the  system. 

For  a  conservative  system  which  has  only  real  eigenvalues  the 
asymptotic  expansion  of  the  Fourier  Inversion  integral  is  obtained 
by  the  technique  of  ray  mathematics.  However,  in  a  non¬ 
conservative  system  the  use  -of  this  technique  combined  with 
analytic  continuation  of  the  eigenvalues  into  the  complex  plane  is 
shown  to  give  incorrect  results  because  the  x/t  mapping  is  not 
unique.  The  correct  expansion  using  the  saddle-point  method 
demonstrates  that  solutions  of  the  dispersion  relation  for  real 
values  of  the  wavenumber  is  not  sufficient.  The  eigenvalues  and 
derivatives  in  the  complex  a-plane  are  required  for  the  correct 
asymptotic  expansion  of  the  initial-value  problem  for  a  non¬ 


conservative  system. 
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APPENDIX  A 


Fourler-Laplace  transform  solution  of  the  initial  value  problem 
In  Section  II  Fourier-Laplace  transforms  are  applied  to  the 
equations  of  motion  to  derive  an  ordinary  differential  equation  in 
the  vertical  perturbation  velocity,  eq.  2.27,  subject  to  the 
boundary  conditions  given  by  eq.  2.19  and  2.20.  Applying  the  Squire 
transformation  to  these  equations  results  in  the  following: 


(s  f  ±aU)  -  a2)  w  -  (s  +  iaU)  UaU") 


+  J(z)  w  =  (s  +  2iay) 


■Jl 

dz 


■  ”o  + 


3w 
_ o 

at 


(A.l) 


=  Q^(a,z;s) 


dw 

dz 


+ 


R 

(s  +  15^)2 


iaU 
_ z 

(s  +  iaU) 


at  z  =  0 


or 


dw  **  ^ 

•^  +  A(oi,s)  w  =  0 


(A.  2) 


and 


w  =  0 


at  z  =  -h 


(A.  3) 
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Dividing  through  the  differential  equation  by  the  coefficient  of  the 
highest  order  terra  gives  the  form 


d^w 

dz2 


+ 


J 

(s  +  iai))2 


±aU" 


(s  +  iai/) 


- 


w 


Q^(a,z,s) 
(s  +  iai/)  2 


(A. 4) 


or 


d^w 

dz2 


+ 


A(J,Z/,s,a) 


Q^(a,z,s) 
(s  +  iat/)2 


Assume  that  (()j(z;s)  and  are  two  independent  solutions 

of  the  homogeneous  differential  equation  such  that 


w(z,s)  =  A<j)^(z,s)  +  B(t)2(z,s) 


(A.  5) 


Making  the  substitution  of  the  solution  in  the  form  A. 5  into  the 
boundary  conditions  A. 2  and  A. 3  leads  to  the  eigenvalue  relation 

F(a,s;J,£/)  =  0  (A.  7) 


for  a  non-trivial  solution  of  the  form  A. 5. 

The  self-adjoint  differential  equation  A. 4  has  the  Green's 
tunction  solution 


w(a,z;s) 


o  G(z,z^)  Q^(a,z,s) 


(s  +  la&)2 


dz 


-h 


(A.  8) 


where  the  Green's  function,  G(z,z^),  is  the  solution  of  the 
differential  equation 


I 
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I 


d^G 

dz2 


+  AG  =  6 (z  + 


(A. 9) 


t 


Let  iJ/j  and  be  two  independent  solutions  of  equation  A. 9  subject 
to  the  following  equations. 

On 

d'I'l 

+  Xi|)j  =0  at  z  =  0  (A.  10) 

On 


ij<2  =  0  at  z  =  -h  (A.  11) 

At  the  point  z  =  the  functions  and  ]p2  ^re  continuous  and 

have  a  jump  in  the  first  derivative  equal  to  unity,  i.e.. 


dij^j 

dz 


z=  -z 


di|»2 

dz 


z=  -z 


(A. 12) 


This  implies  that 


(-z^) 


d<ii. 


1'  o'  dz 


Z"  -7. 


dll', 

-  17 


^  0  (A.  13) 


z=  -z 


where  W  (il'j  ,i|'2)  is  the  Wronskian  of  and  \ti2 
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Applying  the  above  conditions  leads  to  the  following  form  for 
the  Green's  function. 


GCz.Zq)  = 


0  >  z  >  -z 
—  —  o 


(A. 14) 


z  >  z  >  -h 
o  —  — 


The  functions,  can  be  rewritten  in  terms  of  the  linearly 
independent  solutions  of  the  homogeneous  differential  equation, 
A. 4,  <)>j  and  'P^-  The  form  of  the  solution  is 


\li^(z)  =  (|>j(z)(t.2(-h)  -  (t)2(z)(t)j(-h) 

and 


(A. 15) 


d(|)2(o) 

— +X(|>2(o) 


4'2(z) 


d(|>j(o) 

— ^+X<^^(o) 


(A.  16) 


From  these  definitions  of  the  functions,  i); ^ ,  it  follows  that 
the  Wronskian  can  be  written 

W('l'i»'/'2)  “  -F(oiiS)w('<>i  (A. 17) 

and  the  Green's  function  is  of  the  form 
G'(z,z^) 

“  F(5,s)w(^l.<f2) 


(A. 18) 
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Substituting  the  Green's  function  into  equation  A. 8,  gives  the 
solution  to  the  differential  equation  A. 4  in  the  integral  form: 


o 

w(a,2,s)  =  I 
-h 


G'  (z,z^)Q^(a,z,s) 

F(6i,sX«/ ('{’i  .4*2)  (s  +  iay)2 


(A. 19) 


The  function  w(a,z,s)  is  the  Fourier-Laplace  transform  of 
the  solution  to  the  initial  value  problem  which  can  be  written  as 


w(x,z,t)  = 


2iii 


e+i®  o 

1  1 


e-1®  -h 


G' (z,  z^)Q^(a,  z,s) 
F(a,s)w((J)  J  ,<t)2)  (S  +  iaVP 


*  dz  dsda 
o 


(A. 20) 


To  understand  the  information  contained  in  this  equation, 
consider  the  problem  of  a  linear  shear  in  a  homogeneous  fluid, 
i.e.,  y"(z)  =  0  and  J(z)  =  0.  There  are  no  singularities  in  the 
differential  equation  A. 4  and  the  functions  and  ^2  analytic. 
Evaluating  the  integral  over  s  first  in  equation  A. 20  shows  that 
contributions  to  the  integral  are  from  the  poles  of  the  integrand. 
These  contributions  are  in  the  form  of  a  sum  of  exponentials 
arising  from  the  sum  of  the  residues. 

The  poles  of  the  integrand  are  at  the  zeros  of  F(o,s)  and  of 
[s  +  iai/(z^)].  The  function  F(a,s)  =  0  is  the  eigenvalue  relation 
for  the  problem  and  the  contribution  to  the  Integral  from  these 
poles  is  equivalent  to  the  discrete  normal  modes  solution.  The 
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contributions  from  the  zeros  of  [s  +  iai/(z^)]  form  a  continuous 

spectrum  of  modes,  because  the  final  solution  is  also  integrated 

over  z  .  As  Case  (1960)  and  Pedlosky  (1964)  pointed  out,  the 
o 

singular  solutions  to  the  differential  equation  which  arise  from 
the  continuous  spectrum  must  be  added  to  any  normal  modes  solutions 
for  the  initial  value  problem  in  order  to  fit  an  arbitrary  set  of 
initial  conditions. 

If  y"(z)  and  J(z)  are  not  zero  everywhere  in  the  domain  of  the 
problem,  the  situation  is  more  complicated  because  the  differential 
equation  is  singular  at  z^  defined  by 

s  -  iay(z^)  =  0  (A. 21) 

For  simplicity  assume  that  z^  is  a  regular  singular  point.  One 
solution  is  then  of  the  form 

(tij(z)  =  (z  -  z^)*^  Xj  (z)  (A. 22) 

where  Xj (z)  is  analytic. 

The  exponent  r  is  one  of  the  solutions  to  the  indicial 
equation 

r  J(z^)  iaU”  (z^) 

r(r-l)  +  (,  +  [TmsoF  -  jmm  -  “  J  ■  “ 
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If  the  difference  between  the  two  roots,  (r^  -  r2),  is  not  an 
integer,  then  the  other  independent  solution  to  A. 4  is  of  similar 
form: 

^2 

(|>2(z)  =  (z  -  z^)  X2(z) 

However,  if  (r^  -  r2)  is  an  integer  the  form  of  the  solution  is 

(p2(z)  =  (z  -  z^)  lX^(z)ln(z  -  z^)  +  X2(z)]  (A. 25) 

Evaluation  of  the  inverse  Fourier-Laplace  transform  in  eq.  A. 20 
now  includes  the  possibility  of  branch  points  at  the  zeros  of 
(s  +  iaU)  as  well  as  poles. 

This  general  problem  is  discussed  in  a  paper  by  Chimonas  (1979). 
He  again  establishes  the  necessity  of  including  the  continuous 
spectra  of  modes  in 'the  evaluation  of  an  initial  value  problem. 

In  addition,  he  establishes  that  these  modes  may  indeed  lead  to 
instabilities  of  the  system.  They  do  not  appear  as  instabilities 
of  the  vertical  velocity  perturbation,  but  as  horizontal  velocity 
or  density  perturbations  which  grow  algebraically  without  limit. 

In  addition,  he  establishes  the  instability  of  these  solutions  only 
if  the  Richardson  number  of  the  flow  falls  below  0.25  somewhere  in 


the  region  of  the  flow. 


APPENDIX  B 


Asymptotic  Analysis 

The  solution  to  the  initial-value  problem  can  be  evaluated  as 
an  integral  of  the  form: 


w(x,z,t)  =  I  f  (B.l) 

n=l  ■’ 


The  asymptotic  evaluation  of  this  integral  can  be  accomplished  uaing 
the  method  of  steepest  descents.  The  technique  described  below 
follows  the  discussion  in  Jeffreys  and  Jeffreys  (1956). 

Rewrite  the  integral  in  the  form: 

N  * 

w(x,z,t)  =  ^  I  A  w  (a,z)  e*^‘*"^“^da  (B.2) 

n=l  •  -I  " 

where  x/t  -  u»^)  is  a  complex  function  of  the  complex 

variable  a.  The  saddle  point  of  the  function  is  at  a  point 

defined  by 


dd 


0*0. 


0 


or 


x/t  - 


do) 

n 

dd 


-  -  *  0 
0=0^ 


(B.3) 


If  the  second  derivative  of  q  (o)  is  not  zero  at  o^  ,  expand 

n  ” 

q  (o)  in  a  Taylors  series  about  the  saddle  point, 
n 
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Ibe  ijBtecral  ia  S.9  is  of  t3x  form  of  ISatstcm's  lemms.  Tk^ 
first  tern  of  x3ge  asxngttotjLc  ezpaosioxi  of  eq.  1.9  <Bsizt£  Katsm*s 
l«^nniB  is 


wCx.z.tl 


^  W  (L.xl  e^*^^**  e** 

a  • 


^u>  i 

1/2 

■! 

* 

—  'i 

»a  J 

(8-12) 


7^  espasksloB  of  tie  asriBpiotic  analTsis  frcos  tie  two- 
dlneBslonal  case  ^descriled  alorc  to  tie  tlree-diaertsl-CToal  prolleiD 
is  accxxnplisleid  folliowia^  tie  'deriTatiom  Ij  C»ter  (198£|. 

Tie  iatesral  to  he  enralxated  is  sivesi  Ij 


V  (x,J,Z.t) 


■  B 


^i(«»  ♦  Try  -  u^t> 


Cl.  13) 


A;Szizi  ia  tie  asymptotic  lisait  only  tie  most  srcnriac;  vill 

contriltste  to  tie  leadias  terms  of  tie  expaosioo.  He  iate£ral  to 
le  c’valaatcid  is  of  tie  form 

-  I  A^Ci^C«.T.*)  ♦Try  -  ^ 

la  Caster’s  aaalysis  tie  asymptotic  form  of  tie  iate£ral  corner  a  is 


enralaated  first.  For  tlis  iatesratloo.  y  is  treated  as  a  para¬ 
meter.  His  iatepratioa  is  of  tie  same  fora  as  descrilcd  ia  tie 
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two-dimensional  case  and  can  be  evaluated  using  eq.  B.12. 
integral  becomes 


l[a*x/t  +  Yy/t  -  to(a*,Y)]t 

A  W 

\  w(a*,Y,z) 


dY 


d^to 

da2 


(a*.Y) 


1/2 


The 


(B.15) 


where  is  defined  by  eq.  B.3. 

■  Is  ° 


(B.16) 


The  next  step  is  to  expand  the  exponent  in  the  integral  in  eq. 
B.  15  about  the  point  y^-  The  exponent  takes  the  form 


i[a^x/t  +  Y*y/t- u)(a*,Y*)]  t  +  i[y/t--^  (a*,Y*)]  t(Y-Y*) 


if  1  7 

da 

* 

rj  12 
da 

* 

,  (y  - 

1  3y2  "  3ct8Y 

dY 

3a2 

,  dY  J, 

choice  of  y*  to  be 

the 

saddle 

point 

defined  by 

y/t  -  1^  (a*,Y*)  =  0 


(B.17) 


(B.18) 


eliminates  the  linear  terms  in  the  expansion.  Evaluating  the 
integral  as  before  and  making  substitutions  from  above,  the 
asymptotic  form  becomes 


J 
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The  phase  factor  and  constant  have  not  been  carried  along  in 
the  three-dimensional  analysis  because  it  is  the  shape  of  the 
function  that  is  of  primary  interest  and  they  alter  only  the 
amplitude,  not  the  shape. 


b.  Figures  7  and  8 

T_  “  74  cm^/sec^ 


pj  =  1.293  X  10  ^  gm/cm^ 
pj  =  1.02  gm/cm^ 

c.  Figures  9  to  12 
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Figure  5.  Stability  function  (-Q)  vs.  wavenumber 
for  data  set  la.  (i|/  =  angle  between 
mean  velocity  and  wave  normal.) 


Figure  6.  Velocity  ol  stability  boundary  vs. 
wavenumber  for  data  set  la. 


STn#j\.n> 


Figure  7.  Stability  function  (-Q)  vs.  wavenumber 
for  data  set  Ib.  =  angle  between 
mean  velocity  and  wave  normal.) 


Figure  8.  Stability  function  (-Q)  vs.  small  wave- 
number  for  data  set  Ib, 


iNMiNM*  one*  -  nm 


Figure  12.  Asymptotic  expansion  of  Gaussian  pulse 
of  standard  deviation  X  =  0.1. 
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Table  II.  Parameters  for  Figures  13  to  22. 


Set 

'•l 

*^2 

•*3 

h 

U 

s 

U 

m 

A0 

A 

1.0 

2.0 

3.85 

10.0 

0.0 

0.0 

0.0 

0,02 

6 

1.0 

2.0 

3.85 

10.0 

10.0 

10.0 

-.14 

0.0 

C 

1.0 

2.0 

3.85 

10.0 

10.0 

10.0 

-.14 

0.02 

D 

91.0 

92.0 

93.85 

100.0 

10.0 

10.0 

-.14 

0.02 

E 

1.0 

2.0 

3.85 

100.0 

10.0 

10.0 

-.14 

0.02 

F 

1.0 

2.0 

3.85 

'  10.0 

11.0 

10.0 

-.14 

0.02 

G 

1.0 

2.0 

3.85 

10.0 

9.0 

10.0 

-.14 

0.02 

H 

1.0 

7.45 

8.0 

10.0 

10.0 

10,0 

-.14 

0.02 

I 

2.5 

2.53 

3.5 

10.0 

7.11 

7.11 

-.89 

0.02 

J 

2,5 

2.53 

3.5 

10.0 

16.0 

16.0 

-.20 

0.02 

See  Figure  4  for  schematic  diagram  of  profiles  and  parameters 
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Figure  13.  Eigenvalues  for  data  set  A.  No  shear 


0  .5  l.o  US  2.0  2.5  3.0  3.5  4.0  4.5  5.0  5.5  6.0  6.5  7.0  7.5  B.O 

WOVE  NUMBER 

a)  Real  part 


0  .5  1.0  1.5  2.0  2.5  3.0  3.5  4.0  4.5  5.0  5.5  6.0  6.5  7.0  7.5  6.0 

.  ,  .  ,  ,  WAVE  NUMBER 

b)  Imaginary  part 

Figure  14.  Eigenvalues  for  data  set  B.  No  stratification. 
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Figure  15.  Eigenvalues  for  data  set  C.  Combined  profiles. 
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0  .5  1.0  1.5  2.0  2.5  3.0  3.5  4.0  4.5  5.0  5.5  6.0  6.5  7.0  7.5  8.0 

b)  Imaginary  par*  WAVE  NUMBER 

Figure  16.  Eigenvalues  for  data  set  D.  Top 
boundary  far  from  interface. 


IMRGINRRr  OMEGA 


0  .5  1-0  1.5  2.0  2.5  3.0  3.5  4.0  4.5  5.0  5.5  6.0  6.5  7.0  7.5  8.0 

a)  Real  part  WAVE  NUMBER 


0  .5  1.0  1.5  2.0  2.5  3.0  3.5  4.0  4.5  5.0  5.5  6.0  6.5  7.0  7.5  8.0 

b)  Imaginary  part  WAVE  NUMBER 


Figure  17.  Eigenvalues  for  data  set  E.  Bottom 
boundary  far  from  Interface. 


innOH'lR''  OMEGfi 


b)  Itnaglnary  part  WAVE  NUMBER 

Figure  20.  Eigenvalues  for  data  set  H.  Inter¬ 
face  In  lower  half  of  profile. 


Figure  22.  Kigenvnlues  for  dat«  set  J.  High  discharge 


3 

O 


ItlROINfiRY  0ME06 
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b)  Imaginary  part  WflVE—NUfIBEF 


Figure  23.  Eigenv.Tlues  for  Profile  1,  Ap  -  0.002 


IMflGINflRY  OMEGA 


2 


b)  Imaginary  part  WAVE  NUMBER 


Figure  24.  Eigenvalues  for  Profile  1.  Ap  =  0.02 


Figure  25.  F.lgonvalues  for  Profile  1.  Ap  =  0.2. 


IMRCINflRY  OMEGA 
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a)  Real  part  WAVE  NUMBER 


.2 


.1 


0 


-.1 


-.2 

0  .5  1.0  1.5  2.0  2.5  3.0  3.5  4.0  4.5  5.0  5.5  6.0  6.5  7.0  7.5  8.0 

b)  Imaginary  part  WAVE  NUMBER 


Figure  26.  Eigenvalues  for  Profile  2.  Ap  *  0.0 


ItinOINORY  OtiEOfl 


i 


a)  Real  part  WflVE_NUMBE!R 

•2| - 


.1 


0 


1 


0  .5  l.O  1.5  Z.O  ?-.S  3.0  3.5  4..0  4.5  5.0  5.5  6.0  6.5  7.0  7.5  8.0 

b)  Imaginary  part  MOVE  NUMBER 


Figure  27.  r.lgcnval nos  for  Profile  2.  Ap  ■=  0.02 
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0  .5  1.0  1.5  2.0  2.5  3.0  3.5  4.0  4.5  5-0  5.5  6.0  6.5  7.0  7.5  8.0 

b)  Imaginary  part  WAVE  NUMBER 


Figure  28.  Eigenvalues  for  Profile  2.  Ap  =  0.2. 
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b)  Imaginary  part  WAVE!  NUMBER 


Figure  29.  Eigenvalues  for  Profile  3.  Ap  =  0.0. 
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b)  Imaginary  part  WOVE  NUMBER 


Figure  30.  Eigenvalues  for  Profile  3.  Ap  =  0.02. 


Figure  32.  Ei>;cnvaluoR  for  Profile  4.  Ap  =  0.02 


a)  Real  part  WAVE  NUMBER 


Figure  33.  Eigenvalii(?s  for  Profile  b.  Ap  -  0.02. 
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Ap  =  0.02. 


Figure  38.  Time  series.  Initial  exponential  depth  decay 
Profile  3.  a  =  0.5.  (2N  +  6)  modes.  N  -  21 


Contours  of  frequency  in  complex  wave- 
number  plane.  Profile  1.  Ap  ■  0.02. 


Figure  42.  Wavenumber  vs.  x/t  along  line  of  saddle 
points.  Profile  1.  4p  -  0.02. 
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Figure  45.  Asymptotic  expansion.  Initial  flat  spectrum. 

Profile  1,  Ap  =  0.02,  z  =  0.0,  t  =  10.0. 


Figure  47.  Asymptotic  expansion.  Initial  flat  spectrum. 

Profile  1,  Ap  =  0.02,  z  =  -0.2,  t  =  10.0. 
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Figure  48.  Asymptotic  expansion.  Initial  flat  spectrum. 

Profile  1,  Ap  =  0.02,  z  =  -0.5,  t  =  10.0. 
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Figure  49.  Asymptotic  expansion.  Initial  flat  spectrum. 

Profile  1,  Ap  =  0.02,  z  =  0.0,  t  =  100.0. 


••fO  -.10  0  .10  .to  .90  .40  .90  .00  .70  ,90  .10  1 .00  1. 10  l.tO  1.90  |.40  1 .90 

9/T 


Figure  50.  Asymptotic  expansion.  Initial  flat  spectrum. 

Profile  1,  Ao  =  0.02,  z  =  -0.1,  t  =  100.0. 


Figure  52.  Asymptotic  expansion.  Initial  flat  spectrum 
Profile  1,  Ap  =  0.02,  z  =  -0.5,  t  »  100.0. 


Figure  53.  Asymptotic  expansion.  Initial  flat  spectrum.  Second 
mode.  Profile  1,  Ap  =  0.02,  z  =  0.0,  t  =  10.0. 


Figure  54.  Asymptotic  expansion.  Initial  flat  spectrum.  Second 
mode.  Profile  1,  Ap  =  0.02,  z  =  -0.5,  t  =  10.0. 


VCTTIC*.  VOOCITT 


Figure  55.  Asymptotic  expansion.  Initial  flat  spectrum 
mode.  Profile  1,  Ap  =  0.02,  z  =  0.0,  t  =  10( 


Figure  56.  Asymptotic  expansion.  Initial  flat  spectrum 
mode.  Profile  1,  Ap  =■  0.02,  z  =  -0.5,  t  =  1( 


Figure  57.  Asymptotic  expansion.  Initial  Gaussian  spectrum 
Profile  1,  Ap  =  0.02,  z  =  0.0,  t  =  10.0,  X  =  0.1 


Figure  58.  Asymptotic  expansion.  Initial  Gaussian  spectrum 
Profile  1,  Ap  =  0.02,  z  =  0.0,  t  *  100.0,  X  =  0. 
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Figure  59.  Ray  mathematics  expansion.  Initial  flat  spectrum 
Profile  1,  Ap  =  0.02,  z  =  -0.1,  t  =  10.0. 


Figure  60.  Ray  mathematics  expansion.  Initial  flat  spectrum 
Profile  1,  Ap  =  0.02,  z  -  -0.2,  t  =  10.0. 


Figure  61.  Ray  mathematics  expansion.  Initial  flat 
Profile  1,  Ap  =  0.02,  z  =  -0.5,  t  =  10.0. 
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62.  Ray  mathematics  expansion.  Initial  flat  spectrum 
Profile  1,  6p  =  0.02,  z  =  -0.1,  t  =  100.0. 


Figure  63.  Ray  mathematics  expansion.  Initial  flat  spectrum. 
Profile  1,  ftp  =  0.02,  z  »  -0.2,  t  =  100.0. 


64.  Ray  mathematics  expansion.  Initial  flat  spectrum 
Profile  1,  &p  =  0.02,  z  •=  -0.5,  t  =  100.0. 
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Table  IV.  Notation 

a,  b  coefficients 

A,  B  coefficients 

d  non-dimensional  depth  parameter 

f(  )  function 

F(  )  Fourier  transform 

F(  )  dispersion  relation 

g  acceleration  of  gravity 

g(  )  function 

G  body  force  potential 

G(  )  Green’s  function 

h  non-dimensional  depth  of  bottom  boundary 

h^  length  scale 

H  depth  of  bottom  boundary 

1  /T 

j  integer 

J  Richardson  number 

k  integer 

k,  I  horizontal  wavenumber 

L  Laplace  transform 

L  linear  operator 

n  integer 

N  Brunt-VMisMlS  frequency 

p  kinematic  pressure 


w 

I 
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atmospheric  pressure 
q  function 

Q(  )  function 

r  root  of  indicial  equation;  variable 

s  Laplace  transform  variable 

t  time 

u^  velocity  scale 

u,  V  horizontal  perturbation  velocity  components 

U,  V  horizontal  mean  velocity  components 

w  vertical  perturbation  velocity 

W  function 

Wronskian 

X,  y  horizontal  coordinates 

X  function 

z  vertical  coordinate 

a  Y  non-dimensional  horizontal  wavenumbers 

6(  )  Dirac  delta  function 

c  small  parameter 

C  integration  variable 

n  surface  elevation 

6  non-dimensional  density  perturbation;  angle 

6  non-dimensional  mean  density 

X  standard  deviation;  parameter;  function 


A 


function 
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p  density 

density  scale 
o  frequency 

T  coefficient  of  surface  tension 

^  function 

t|)  angle  between  wavefront  normal  and  mean  velocity; 

function 

“  non-dimensional  frequency 

Symbols: 

(  )  initial  value;  Initial  condition 

o 

(  )^  saddle  point  value;  singular  point 
(  )g  particular  solution  a 

(  ),  particular  solution  b 

D 

(  profile  I 

(  )j,j.  profile  II 

(  profile  III 

(  variable  in  j-th  layer 

(  )  n-th  mode 

n 

(_)  dimensional  variable 

t 

(_)  mean  plus  perturbation  variable 

(  )  mean  value 

I 

(  )  derivative  with  respect  to  z 

(  )"  second  derivative  with  respect  to  z 

(  )  complex  conjugate 

(  )  Fourier  transformed  variable;  Squire  transformed  variable 
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